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ABSTRACT 

Finite volume methods are a class of discretization schemes that have proven highly successful in 
approximating the solution of a wide variety of conservation law systems. They are extensively used 
in fluid mechanics, porous media flow, meteorology, electromagnetics, models of biological processes, 
semi-conductor device simulation and many other engineering areas governed by conservative systems 
that can be written in integral control volume form. 

This article reviews elements of the foundation and analysis of modern finite volume methods. 
The primary advantages of these methods are numerical robustness through the obtention of discrete 
maximum (minimum) principles, applicability on very general unstructured meshes, and the intrinsic 
local conservation properties of the resulting schemes. Throughout this article, specific attention is 
given to scalar nonlinear hyperbolic conservation laws and the development of high order accurate 
schemes for discretizing them. A key tool in the design and analysis of finite volume schemes suitable for 
non-oscillatory discontinuity capturing is discrete maximum principle analysis. A number of building 
blocks used in the development of numerical schemes possessing local discrete maximum principles are 
reviewed in one and several space dimensions, e.g. monotone fluxes, E-fluxes, TVD discretization, non- 
oscillatory reconstruction, slope limiters, positive coefficient schemes, etc. When available, theoretical 
results concerning a priori and a posteriori error estimates are given. Further advanced topics are then 
considered such as high order time integration, discretization of diffusion terms and the extension to 
systems of nonlinear conservation laws. 
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maximum principles, higher order schemes 
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1. Introduction: Scalar nonlinear conservation laws 

Many problems from physics, chemistry, biology, mechanics, and gas dynamics lead to the 
study of nonlinear hyperbolic conservation laws. As a prototype conservation law, consider the 
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Cauchy initial value problem 

d t u 4- V • f(u) = 0 in R d x E + , (la) 

u(x, 0) = u 0 (rr) in E d . (lb) 

Here u(x,t) : l' 1 x 1 + — > E denotes the dependent solution variable, f(u) € 6 fl (M) denotes 
the flux function, and uq(x) : E d -> E the initial data. 

The function u is a classical solution of the scalar initial value problem if u e C 1 (E d x E + ) 
satisfies (la-lb) pointwise. An essential feature of nonlinear conservation laws is that, in 
general, gradients of u blow up in finite time, even when the initial data uo is arbitrarily 
smooth. Beyond some critical time to classical solutions of (la-lb) do not exist. This behavior 
will be demonstrated using the method of characteristics. By introducing the notion of weak 
solutions of (la-lb) together with an entropy condition, it then becomes possible to define 
a class of solutions where existence and uniqueness is guaranteed for times greater than to- 
These are precisely the solutions that are numerically sought in the finite volume method. 

LI. Characteristics of scalar conservation laws 

Let u be a classical solution of (la) subject to initial data (lb). Further, define the vector 

o(u) = /'(«) = ,fi(u)) T . 

A characteristic T y is a curve (x(t),t) such that 

x f (t) = a(u(x(t),t)) for t > 0 , 

x(0) = y . 

Since u is assumed to be a classical solution, it is readily verified that 

^-u(x(t), t) = d t u + x'(t)Vu = d t u + a(u) Vu = d t u -h V ■ f(u) = 0 . 
at 

Therefore, u is constant along a characteristic curve and F y is a straight line since 
x f (t) = a(u(x(t),t)) = a(u(x(0),0)) = a(u(y, 0)) = a(u 0 (y)) = const . 

In particular x(t) is given by 

x(t) = y + ta(u 0 (y)) • (2) 

This important property may be used to construct classical solutions. If z is fixed and y 
determined as a solution of (2), then 

u(x,t) = u 0 (y) . 

This procedure is the basis of the so-called method of characteristics. On the other hand, 
this construction shows that the intersection of any two straight characteristic lines leads to a 
contradiction in the definition of u(x,t). Thus, classical solutions can only exist up to the first 
time at which any two characteristics intersect. 
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1.2. Weak solutions 


Since, in general, classical solutions only exist for a finite time t 0 , it is necessary to introduce 
the notion of weak solutions that are well-defined for times t > to. 

Definition 1.1 (Weak solution) Let u 0 € L°°( R d ). Then , u is a weak solution of (la-lb) 
ifu€ L°°(R d x R + ) and (la-lb) hold in the distributional sense , i.e. 


J j (ud t <f> + /(u)-V^) 

R d R+ 


dt dx + / uo<f>(x , 0) dx = 0 
R d 


for all (j> e Cl(R d x R+) . 


(3) 


Note that classical solutions are weak solutions and weak solutions that lie in C 1 (R d x R + ) 
satisfy (la-lb) in the classical sense. 

It can be shown (see Kruzkov, 1970; Oleinik, 1963) that there always exists at least one 
weak solution to (la-lb) if the flux function / is at least Lipschitz continuous. Nevertheless, 
the class of weak solutions is too large to ensure uniqueness of solutions. An important class of 
solutions are piecewise classical solutions with discontinuities separating the smooth regions. 
The following lemma gives a necessary and sufficient condition imposed on these discontinuities 
such that the solution is a weak solution (see for example Godlewski and Raviart, 1991; Kroner, 
1997). Later a simple example is given where infinitely many weak solutions exist. 

Lemma 1.2 (Rankine-Hugoniot jump condition) Assume that R d x R + is separated by 
a smooth hypersurface S into two parts Qi and Q r - Furthermore , assume u is a C l -function 
on fi* and respectively. Then , u is a weak solution of (la-lb) if and only if the following 
two conditions hold: 

i) u is a classical solution in £li and fl r , respectively. 

ii) u satisfies the Rankine- Hug oniot jump condition , i.e. 


[u]s = [f(u)] * v on S . (4) 

Here , ( i /, — s) T denotes a unit normal vector for the hypersurface S and [w] denotes the jump 
in w across the hypersurface S . 


In one space dimension, it may be assumed that S is parameterized by (a(t),t) such that 
5 — cr ; (t) and v = 1. The Rankine-Hugoniot jump condition then reduces to 


s = 


[/(«)] 

M 


on 5 . 


(5) 


Example 1.3 (Non- uniqueness of weak solutions) Consider the one- dimensional Burg- 
ers 7 equation , f(u) = u 2 / 2, with Riemann data: Uo(x) = ui for x < 0 and uo(x) = u r for 
x > 0. Then, for any a > max(ui , — u T ) a function u given by 


u(x,t) = < 


Hi, 

X < Sit 

-a, 

Sit < x < 0 

a, 

0 < X < S2t 

u r , 

s 2 t < X 

S2 = 

(a -h u r )j 2. 


( 6 ) 


piecewise constant and satisfies the Rankine-Hugoniot jump condition. This elucidates a one- 
parameter family of weak solutions. In fact, there is also a classical solution whenever ui < u r . 
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In this case , the characteristics do not intersect and the method of characteristics yields the 
classical solution 

{ ui , x < U\t 

x/t, uit < x < u r t . (7) 

u r , u r t < x 

This solution is the unique classical solution but not the unique weak solution. Consequently, 
additional conditions must be introduced in order to single out one solution within the class 
of weak solutions. These additional conditions give rise to the notion of a unique entropy weak 
solution. 

1.3. Entropy weak solutions and vanishing viscosity 

In order to introduce the notion of entropy weak solutions, it is useful to first demonstrate 
that there is a class of additional conservation laws for any classical solution of (la). Let u be 
a classical solution and tj : R R a smooth function. Multiplying (la) by 7]'(u ), one obtains 

0 — r} f (u)d t u 4* r) f (u)V • f(u) = d t r)(u) -f V • F(u) (8) 

where F is any primitive of 7 //'. This reveals that for a classical solution u, the quantity r)(u ), 
henceforth called an entropy function, is a conserved quantity. 

Definition 1.4 (Entropy - entropy flux pair) Let r] : E -4 M be a smooth convex function 
and F : R R a smooth function such that 

F' = if (9) 

in (8). Then (77, F) is called an entropy - entropy flux pair or more simply an entropy pair for 
the equation (la). 

Note 1.5 (Kruzkov entropies) The family of smooth convex entropies r} may be 
equivalently replaced by the non-smooth family of so-called Kruzkov entropies , i.e. r] K (u) = 
ju — k\ for all t c € R ( see Kroner , 1997). 

Unfortunately, the relation (8) can not be fulfilled for weak solutions in general, as it would lead 
to additional jump conditions which would contradict the Rankine-Hugoniot jump condition 
lemma. Rather, a weak solution may satisfy the relation (8) in the distributional sense with 
inequality. To see that this concept of entropy effectively selects a unique, physically relevant 
solution among all weak solutions, consider the viscosity perturbed equation 

d t u € 4 V • f(u e ) = e A u e (10) 

with e > 0. For this parabolic problem, it may be assumed that a unique smooth solution u e 
exists. Multiplying by rf and rearranging terms yields the additional equation 

d t r)(u e ) + V • F(u t ) = eA r?(u e ) - €7?"(it e )|Vu| 2 . 

Furthermore, since rj is assumed convex (rf* > 0), the following inequality is obtained 

dtV(Ue) 4* V * F(u € ) < eAr}(u e ) . 

Taking the limit e -4 0 establishes (see Malek, Necas, Rokyta and Ruzicka, 1996) that u € 
converges towards some u a.e. in x R + where u is a weak solution of (la-lb) and satisfies 
the entropy condition 

d t rj(u) + V-F(u)< 0 (11) 
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in the sense of distributions on R d xl + . 

By this procedure, a unique weak solution has been identified as the limit of the 
approximating sequence u t . The obtained solution u is called the vanishing viscosity weak 
solution of (la-lb). Motivated by the entropy inequality (11) of the vanishing viscosity solution, 
it is now possible to introduce the notion of entropy weak solutions. This notion is weak enough 
for the existence and strong enough for the uniqueness of solutions to (la-lb). 

Definition 1.6 (Entropy weak solution) Let u be a weak solution of (la- lb). Then , u is 
called an entropy weak solution if u satisfies for all entropy pairs ( 77 , F) 

(r)(u)d t <t> 4- F(u) • V<f>) dtdx + 

R d R + R d 

From the vanishing viscosity method, it is known that entropy weak solutions exist. The 
following L 1 contraction principle guarantees that entropy solutions are uniquely defined (see 
Kruzkov, 1970). 

Theorem 1.7 (L 1 contraction principle) Let u and v be two entropy weak solutions of 
(la- lb) with respect to initial data uq and v$. Then , the following L 1 contraction principle 
holds 

M-,f) -u(-,f)IUl(R<») < ||« 0 -Vo||il(R-). (13) 

for almost every t > 0. 

This principle demonstrates a continuous dependence of the solution on the initial data and 
consequently the uniqueness of entropy weak solutions. Finally, note that an analog of the 
Rankine-Hugoniot condition exists (with inequality) in terms of the entropy pair for all entropy 
weak solutions 

[r)(u)] s > [F(u)] * v on S . (14) 

1.4- Measure-valued or entropy process solutions 

The numerical analysis of conservation laws requires an even weaker formulation of solutions to 
(la-lb). For instance, the convergence analysis of finite volume schemes makes it necessary to 
introduce so called measure-valued or entropy process solutions (see DiPerna, 1985; Eymard, 
Galluoet and Her bin, 2000). 

Definition 1.8 (Entropy process solution) A function p(x ,t,a) E L CX) (R d xl + x (0,1)) 
is called an entropy process solution of (la- lb) if u satisfies for all entropy pairs ( 77 , F) 

1 

T){p)d t (<j> 4- F(fj) • V0) da dt dx 4- J r](uo)<p(x, 0) dx > 0 for all <fi £ Cq (R d x R + , E + ) . 

RdR + 0 Rd 

The most important property of such entropy process solutions is the following uniqueness and 
regularity result (see Eymard, Galluoet and Herbin, 2000 [Theorem 6.3]). 

Theorem 1.9 (Uniqueness of entropy process solutions) Let uq £ L°°( R d ) and f E 
G 1 (M). The entropy process solution p of problem (la-lb) is unique . Moreover , there exists a 
function u E L°°(R d xl + ) such that u(x,t) = p(x,t,a) a.e. for (x,t,a) E R d x 1+ x (0,1) 
and u is the unique entropy weak solution of (la- lb). 




f rj(uo)<j){x, 0) dx > 0 for all <f> 6 C% (R d xl + ,l + ). (12) 
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2. Finite volume (FV) methods for nonlinear conservation laws 
In the finite volume method, the computational domain, Q C is first tessellated into a 
collection of non overlapping control volumes that completely cover the domain. Notationally, 
let T denote a tessellation of the domain fi with control volumes T £ T such that UtgtT = fh 
Let hr denote a length scale associated with each control volume T , e.g. hr = diam(T). For 
two distinct control volumes T; and Tj in T, the intersection is either an oriented edge (2-D) 
or face (3-D) with oriented normal or else a set of measure at most d— 2. In each control 
volume, an integral conservation law statement is then imposed. 


Definition 2.1 (Integral conservation law) An integral conservation law asserts that the 
rate of change of the total amount of a substance with density u in a fixed control volume T is 
equal to the flux f of the substance through the boundary dT 


d_ 

dt 


L 


udx + 



• dv = 0 . 


(15) 


This integral conservation law statement is readily obtained upon spatial integration of the 
divergence equation (la) in the region T and application of the divergence theorem. The 
choice of control volume tessellation is flexible in the finite volume method. For example, Fig. 



a. Cell-centered 


b. Vertex-centered 


Figure 1. Control volume variants used in the finite volume method: 
(a) cell-centered and (b) vertex-centered control volume tessellation. 


1 depicts a 2-D triangle complex and two typical control volume tessellations (among many 
others) used in the finite volume method. In the cell-centered finite volume method shown in 
Fig. la, the triangles themselves serve as control volumes with solution unknowns stored on 
a per triangle basis. In the vertex- centered finite volume method shown in Fig. lb, control 
volumes are formed as a geometric dual to the triangle complex and solution unknowns stored 
on a per triangulation vertex basis. 
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2.1. Godunov finite volume discretizations 

Fundamental to finite volume methods is the introduction of the control volume cell average 
for each Tj G T 



For stationary meshes, the finite volume method can be interpreted as producing an evolution 
equation for cell averages 

FtL l ndx = lTil Ft Ut ■ (I7) 

Godunov, 1959 pursued this interpretation in the discretization of the gas dynamic equations 
by assuming piecewise constant solution representations in each control volume with value 
equal to the cell average. However, the use of piecewise constant representations renders the 
numerical solution multivalued at control volume interfaces thereby making the calculation of 
a single solution flux at these interfaces ambiguous. The second aspect of Godunov’s scheme 
and subsequent variants was the idea of supplanting the true flux at interfaces by a numerical 
flux function , g(u,v) a Lipschitz continuous function of the two interface states 

u and v. A single unique numerical flux was then calculated from an exact or approximate 
local solution of the Riemann problem in gas dynamics posed at these interfaces. Figure 2 
depicts a representative 1-D solution profile in Godunov’s method. For a given control volume 
Tj — [xj_i/2,Xj+i/ 2 \, Riemann problems are solved at each interface Xj± i/ 2 . For example, at 
the interface Xj +1 f 2 the Riemann problem counterpart of (la-lb) 

d T W j+1 / 2 {Z,T) + 0c/(tw,- + i /2 (f,T)) = 0 in E x E + 
for Wj+ 1/2 (£,t) € K with initial data 

m;jj 

is solved either exactly or approximately. From this local solution, a single unique numerical 
flux at Xj+ 1/2 is computed from g(uj,Uj+ 1 ) = 2 (0,E + )). This construction utilizes the 

fact that the solution of the Riemann problem at £ = 0 is a constant for all time r > 0. 

u(x,t) 


X j-3/2 Xj-1/2 X j+l/2 Xj +3 /2 

Figure 2. 1 -D control volume, Tj = [£7-1/2, £7+1/2], depicting Godunov’s interface 
Riemann problems, wj± 1/2 from piecewise constant interface states. 
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In higher space dimensions, the flux integral appearing in (15) is similarly approximated by 



■dvK, ^2 9jk{uj,u k ) 

veneer* 


(18) 


where the numerical flux is assumed to satisfy the properties: 


• (Conservation) This property insures that fluxes from adjacent control volumes sharing 
a mutual interface exactly cancel when summed. This is achieved if the numerical flux 
satisfies the identity 

9jk ( u , v ) = -g k j (v,u) . (19a) 

• (Consistency) Consistency is obtained if the numerical flux with identical state arguments 
reduces to the true flux of that same state, i.e. 

g jk (u,u)=f f (u) • dv . (19b) 

Je jh 

Combining (17) and (18) yields perhaps the simplest finite volume scheme in semi-discrete 
form. Let V£ denote the space of piecewise constants, i.e. 

V° = {« | v\t € x(T), VT € T) (20) 

with x(T) a characteristic function in the control volume T. 

Definition 2.2 (Semi-discrete finite volume method) The semi-discrete finite volume 
approximation of (la- lb) utilizing continuous in time solution representation, t 6 [0, r], and 
piecewise constant solution representation in space, Uh(t ) € V£, such that 


with initial data 


Uj(t ) = p— f U h {x,t ) dx 

Mil Jt, 

u j (0) = u 0 {x) dx 


and numerical flux function gjk(uj, Uk) is given by the following system of ordinary differential 
equations 

^ E 3jk (Uj, U k ) = 0 , VT ) € r . (21) 


Ve jh edTj 


This system of ordinary differential equations can be marched forward using a variety of 
explicit and implicit time integration methods. In Sect. 4.1, time integration schemes that 
preserve properties of the spatial discretization are considered in more detail. Let uf denote a 
numerical approximation of the cell average solution in the control volume Tj at time t n = nAL 
A particularly simple time integration method is the forward Euler scheme 


d 

— Uj 

dt 3 


u T l ~ u l 

At 


thus producing the fully-discrete finite volume form. 
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Definition 2.3 (Fully-discrete finite volume method) The fully- discrete finite volume 
approximation of (la- lb) for the time slab interval [t n ,t n A- At] utilizing the piecewise constant 
solution representation in space, u% € V£, such that 


u ‘ = \k 


with initial data 


= j^i J TI “" (i ) ix 

and numerical flux function gjkiu'j , u£) is given by the following fully-discrete system 


„n+i 


A t 


= T 9 jk (u],u n k ) , VT,- e T . 


( 22 ) 


Ji Vejfc £dTj 


In subsequent sections, properties of the semi-discrete scheme (21) and fully-discrete scheme 
(22) will be examined in more detail. 


2.1.1. Monotone schemes. Unfortunately, the numerical flux conditions (19a) and (19b) are 
insufficient to guarantee convergence to entropy satisfying weak solutions (12) and additional 
numerical flux restrictions are necessary. Two classes of numerical fluxes that guarantee such 
convergence for piecewise constant numerical solution data are monotone fluxes and E-fluxes. 
Specifically, Harten, Hyman and Lax, 1976 provide the following result concerning convergence 
of the fully-discrete one-dimensional scheme to weak solutions which was later generalized to 
(22) and irregular grids by Cockburn, Coquel and Lefloch, 1994. 


Theorem 2.4 (Monotone schemes and weak solutions) Consider a 1-D finite volume 
discretization of (la- lb) with 2k A - 1 stencil on a uniformly spaced mesh in both time and space 
with corresponding mesh spacing parameters At and Ax 

Uj = Hj {Uj+k , - • • , Uj , . . . , Uj—k) 

= U 1 - ^ -{dj+ 1/2 - 9j- 1 / 2 ) (23) 

and consistent numerical flux of the form 


9j+ 1/2 ~ d( u j+k ? • • • 3 u j+li Uj, ... , Uj-k- j-l) 


that is monotone in the sense 


dHj 

duj+i 


> 0 


< k 


(24) 


Then as At and Ax tend to zero with At/ Ax = constant, u ? converges boundedly almost 
everywhere to u(x,t), an entropy satisfying weak solution of (la-lb). 


The monotonicity condition (24) motivates the introduction of Lipschitz continuous monotone 
fluxes satisfying 


> 0 if i=j 

OUi 

%!£< 0if if. 


(25a) 

(25b) 
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together with a CFL (Courant-Friedrichs-Levy) like condition 

x At / dfr+i/2 ^-i/2 \ >0 
Ax \ duj duj ) ” 

so that (24) is satisfied. Some examples of monotone fluxes for (la) include 


• (Godunov flux) 


• (Lax-Friedrichs flux) 


min t{u) 

__ J ue[u j ,u j + 1 ] 

| max f(u) 

{ue[uj,u 3+ 1] 


if uj < Uj+x 
if Uj > Uj+i 


9jh/2 - 9 + /( w i+l)) “ 9 SU P I /»! ( u i+l “ «i) • 

1 1 «e[ui,-ui+i] 


(26) 


(27) 


2.L2. E-flux schemes. A more general class of numerical fluxes was introduced and analyzed 
by Osher, 1984 that still guarantees convergence to weak entropy solutions when used in (22) 
or (23). These fluxes are called E-fluxes, gj+ 1/2 = g E {uj+k , . • • ,Uj+i,Uj, * * • ,%-fc+i), due to 
the relationship to Olienick’s well-known E-condition which characterizes entropy satisfying 
discontinuities. E-fluxes satisfy the inequality 

gf +1/2 - /(«) 0 Vu € [ Uj>u ] . (28) 

Uj + 1 - Uj 

E-fluxes can be characterized by their relationship to Godunov’s flux. Specifically, E-fluxes are 
precisely those fluxes such that 

df+i /2 — $+ 1/2 u j + 1 ^ u j (29a) 

df+1/2 > 9?+l/2 if 1 > U 3 ■ ( 29b ) 

Viewed another way, note that any numerical flux can be written in the form 

9j+ 1/2 = ^ (/(«;) + f{uj+i)) - ^Q(uj + 1 - uj) (30) 

where Q(-) denotes a viscosity for the scheme. When written in this form, E-fluxes are those 
fluxes that contribute at least as much viscosity as Godunov’s flux, i.e. 

Qf+ 1/2 ^ Qj+ 1/2 • ( 31 ) 

The most prominent E-flux is the Enquist-Osher flux 

1 1 f U 3 + 1 

gf ?!, 2 = I (/(«;) + /(%•+!)) ~ 2 l/'( s )l ds ’ ( 32 ) 

although other fluxes such as certain forms of Roe’s flux with entropy fix fall into this category. 
From (29a-29b), the monotone fluxes of Godunov gf +1 / 2 an <i Lax-Friedrichs 9 ^+^/ 2 are also 
E-fluxes. 
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2.2. Stability , convergence and error estimates 

Several stability results are provided here that originate from discrete maximum principle 
analysis and are straightforwardly stated in multi-dimensions and on general unstructured 
meshes. In presenting results concerning convergence and error estimates, a notable difference 
arises between one and several space dimensions. This is due to the lack of a BV bound on 
the approximate solution in multi-dimensions. Thus, before considering convergence and error 
estimates for finite volume methods, stability results are presented first together with some a 
priori bounds on the approximate solution. 


2.2.1. Discrete maximum principles and stability. A compelling motivation for the use of 
monotone and E-fluxes in the finite volume schemes (21) and (22) is the obtention of discrete 
maximum principles in the resulting numerical solution of nonlinear conservation laws (la). A 
standard analysis technique is to first construct local discrete maximum principles which can 
then be applied successively to obtain global maximum principles and stability results. 

The first result concerns the boundedness of local extrema in time for semi-discrete finite 
volume schemes that can be written in nonnegative coefficient form. 

Theorem 2.5 (LED Property) The semi-discrete scheme for each Tj G T 

~TT = TjTT Cjk{uh)(uk — Uj), (33) 

1 Ve jk €dTj 

is Local Extremum Diminishing (LED), i.e. local maxima are non-increasing and local minima 
are nondecreasing , if 

Cjk(u>h) > 0 j Vej* G d Tj . (34) 


Rewriting the semi-discrete finite volume scheme (21) in terms of monotone fluxes or E-fluxes 

9jk {Uj , Uk) — f (Uj) • Vjk / —xl A 
u k - uj k 3 

= -]4t ~ u i) ( 35 ) 

1 - 7 ’ 1 Ve jk edTi k 

for appropriately chosen Uj k e [uj, u k ] together with the monotone flux conditions (25a-25b) 
or the E-flux condition (28) reveals that monotone flux and E-flux finite volume schemes 
(21) are LED. In order to obtain local space-time maximum principle results for the fully- 
discrete discretization (22) requires the introduction of an additional CFL-like condition for 
non-negativity of coefficients in space-time. 

Theorem 2.6 (Local space-time discrete maximum principle) The fully- discrete scheme 
for the time slab increment [i n ,£ n+1 ] and each Tj G T 


du j 1 \ 

Ik = ”i t7\ ^ 


\Tj\ 


Ve^fc EdT. 


u]+'=u? + — C jk (K)(u n k -u») 

1 3 ' VejkZdT; 


exhibits a local space-time discrete maximum principle 


min 

Ve^eST,- 


(«*. 


<) < < +1 


< max ( u n k ,uJ) 
k € oTj 


(36) 


(37) 
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if 

c jk K)>0, Ve iJfe e dTj 

and satisfies the CFL-like condition 

i-#r E CM<)>0. 

1 Ve jh £.dTj 


(38) 

(39) 


Again noting that the flux terms in the fully-discrete finite volume scheme (22) can be written 
in the form (35) reveals that the monotone flux conditions (25a-25b) or the E-flux condition 
(28) together with a local CFL-like condition obtained from (39) imply a local space-time 
discrete maximum principle. By successive application of Theorem 2.6, a global Testability 
bound is obtained for the scalar initial value problem (la-lb) in terms of initial data uo(x). 

Theorem 2.7 (jL°°-stability) Assume a fully-discrete finite volume scheme (22) for the 
scalar initial value problem (la- lb) utilizing monotone fluxes or E- fluxes that satisfy a local 
CFL-like condition as given in Theorem 2.6 for each time slab increment [t n ,£ n-h1 ]. Under 
these conditions , the finite volume scheme is L°° -stable and the following estimate holds: 

inf uo(:r) <u]< sup uq{x ), for all ( Tj,t n ) € T x [Q,r]. (40) 

xeR* xeR d 

Consider now steady-state solutions, u n+1 = u n = u*, using monotone flux or E-flux schemes in 
the fully-discrete finite volume scheme (22). At steady state, noil-negativity of the coefficients 
C(uh ) in (36) implies a discrete maximum principle. 

Theorem 2.8 (Local discrete maximum principle in space) The fully-discrete scheme 
(36) exhibits a local discrete maximum principle at steady state, u for each Tj € T 

min ul < u*j < max ut (41) 

\/e jk edTj * 


if 

Cjk(uh) > 0 , Ve^-fc € dTj . 

Once again by virtue of (25a-25b) and (28), the conditions for a local discrete maximum 
principle at steady state are fulfilled by monotone flux and E-flux finite volume schemes (22). 
Global maximum principles for characteristic boundary valued problems are readily obtained 
by successive application of the local maximum principle result. 

The local maximum principles given in (37) and (41) preclude the introduction of spurious 
extrema and 0(1) Gibbs-like oscillations that occur near solution discontinuities computed 
using many numerical methods (even in the presence of grid refinement). For this reason, 
discrete maximum principles of this type are a highly sought after design principle in the 
development of numerical schemes for nonlinear conservation laws. 


2.2.2. Convergence results. The Testability bound (40) is an essential ingredient in the 
proof of convergence of the fully-discrete finite volume scheme (22). This bound permits the 
subtraction of a subsequence that converges against some limit in the T°° weak-* sense. The 
primary task that then remains is to identify this limit with the unique solution of the problem. 
So although Testability is enough to ascertain convergence of the scheme, stronger estimates 
are needed in order to derive convergence rates. 
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Let BV denote the space of functions with bounded variation, i.e. 

BV = {g € L 1 ^) | \g\sv < 00 } with \g\sv = sup [ g V ■ <p dx . 

JR d 

From the theory of scalar conservation laws, it is known that, provided the initial data is in 
BV, the solution remains in BV for all times. Therefore, it is desirable to have an analog of 
this property for the approximate solution as well. Unfortunately, up to now, such an result is 
only rigorously proved in the one-dimensional case or in the case of tensor product cartesian 
meshes in multiple space dimensions. In the general multi-dimensional case, the approximate 
solution can only be shown to fulfill some weaker estimate which is thus called a weak BV 
estimate (see Vila, 1994; Cockburn, Coquel and Lefloch, 1994; Eymard, Gallouet, Ghilani and 
Herbin, 1998). 

Theorem 2.9 (Weak BV estimate) Let T be a regular triangulation , and let J be a 
uniform partition of [0,r], e.g. A t n = At. Assume that there exists some a > 0 such that 
ah 2 < |T*;j, a < h. For the time step A t n , assume the following CFL-like condition for 
a given £ € (0, 1) 

A." < (V-O- 3 * 

Lg 

where L g is the Lipschitz constant of the numerical flux function . Furthermore , let uq 6 
L°°(W i ) fl BV(R d ) f! L 2 (R d ). Then, the numerical solution of the fully- discrete discretization 
(22) fulfills the following estimate 

K-uriQ,-, («".«") < KVT\B R+h (0)\Vh , (42) 

n jl 

where K only depends on a, L g , £ and the initial function uq. In this formula Qji is defined 
as 

„ , X _ 2 9ji(u, v ) - gji(u, u ) - gji(v, v) 

Qjl{UiV) = 

and Xji denotes the discrete cutoff function on Br(0 ) C R d , i.e. 

/ 1, 2/ (Tj U Ti) D Br(0) 0 

Xjl “ \ 0, else 

Note that in the case of a strong BV estimate, the right-hand side of (42) would be 0(h) 
instead of 0(Vh). 

Another important property of monotone finite volume schemes is that they preserve the 
I/ 1 -contraction property (see Theorem 1.7). 

Theorem 2.10 (I^-contraction property and Lipschitz estimate in time) Letuh.Vh € 
V^ 0 be the approximate monotone finite volume solutions corresponding to initial data uq, Vo 
assuming that the CFL-like condition for stability has been fulfilled. Then the following discrete 
L 1 -contraction property holds 

IK(*,t + r) + r)|| jL i (R <i ) < |K(*,t) -tffc(-,0lli,i(R<9 • 

Furthermore , a discrete Lipschitz estimate in time is obtained 

E l T ;IK +1 - ^ L ^ tn E E l e i'lK - u i\ • 

j J l 
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The principle ingredients of the convergence theory for scalar nonlinear conservation laws 
are compactness of the family of approximate solutions and the passage to the limit within 
the entropy inequality (12). In dealing with nonlinear equations, strong compactness is needed 
in order to pass to the limit in (12). In one space dimension, due to the BV estimate and 
the selection principle of Helly, strong compactness is ensured and the passage to the limit is 
summarized in the well known Lax- Wendroff theorem (see Lax and Wendroff, 1960). 

Theorem 2.11 ( Lax- Wendroff theorem) Let (u m ) m€ N be a sequence of discrete solutions 
defined by the finite volume scheme in one space dimension with respect to initial data u$. 
Assume that (u m )m€N is uniformly bounded with respect to m in L°° and u m converges almost 
everywhere in R x R + against some function u. Then u is the uniquely defined entropy weak 
solution of (la- lb). 

With the lack of a BV estimate for the approximate solution in multiple space dimensions, 
one cannot expect a passage to the limit of the nonlinear terms in the entropy inequality in 
the classical sense, i.e. the limit of u m will not in general be a weak solution. Nevertheless, 
the weak compactness obtained by the Z°°-estimate is enough to obtain a measure-valued or 
entropy process solution in the limit. 

The key theorem for this convergence result is the following compactness theorem of Tartar 
(see Tartar, 1983; Eymard, Galluoet and Herbin, 2000). 

Theorem 2.12 (Tartar’s Theorem) Let (u m ) mG N be a family of bounded functions in 
L°°(E n ). Then, there exists a subsequence (u m ) m€ N, and a function u E L°°(M n x (0,1)) 
such that for all functions g G C7(1E) the weak-* limit of g(u m ) exists and 

lim f g(u m (x))<j>(x)dx = f [ g(u(x,a))(f)(x)dx da, for all <f G L 1 (R n ) . (43) 

m^°° J Rn J 0 J Rn 

In order to prove the convergence of a finite volume method, it now remains to show that 
the residual of the entropy inequality (12) for the approximate solution Uh tends to zero if h 
and At tend to zero. Before presenting this estimate for the finite volume approximation, a 
general convergence theorem is given which can be viewed as a generalization of the classical 
Lax- Wendroff result (see Eymard, Galluoet and Herbin, 2000). 

Theorem 2.13 (Sufficient condition for convergence) Let uo G L°°(R d ) and f E 
C fl (E). Further , let (u m ) m€ N be any family of uniformly bounded functions in L°°(R d x ) 
that satisfies the following estimate for the residual of the entropy inequality using the class of 
Kruzkov entropies rj K ( see Note 1.5). 

/ / {^ K [um)dt4> + F VK (u m ) * V<t>) dtdx + / 7] K (uo)(p{x,0)dx > -R(K,u m ,</>) (44) 

jR d JR+ JU. d 

for all k GR and <f> G Cq (® d x M + , ) where the residual R(k , u m , (j>) tends to zero for m — > oo 

uniformly in k. Then , u m converges strongly to the unique entropy weak solution of (la-lb) in 
Lf oc (R d x R+) for all p G [1, oo) . 

Theorem 2.14 (Estimate on the residual of the entropy inequality) Let (u m ) m€ N be 
a sequence of monotone finite volume approximations satisfying a local CFL-like condition 
as given in (39) such that h,At tend to zero for m oo. Then , there exist measures 
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ix m £ M(R d xl + ) and v m £ A4(R^) such that the residual R(K,u m , (f) of the entropy inequality 
is estimated by 

R{n,u m ,<f>)< f [ {\d t <j>{x y t)\ + \V<j>(x,t)\)dp, m (x,t) + [ 0)du m (x) 

J R d J R+ J R d 

for all k £ R and <fi £ Co(M d x R + ,R + ). The measures p m and v m satisfy the following 
properties: 

L For all compact subsets ft CC R d x R+ , lim™-^ /i m (ft) = 0. 

2 . For all g £ Co(M d ) the measure u m is given by { v m ,g ) — J Rd g(x)\u 0 {x) - u m (x,0)\dx. 

These theorems are sufficient for establishing convergence of monotone finite volume schemes. 

Corollary 2.15 (Convergence theorem) Let (u m ) meN be a sequence of monotone finite 
volume approximations satisfying the assumptions of Theorem 2.14 • Then , u m converges 
strongly to the unique entropy weak solution of (la- lb) in Lf oc (R d x R + ) for all p £ [l,oo). 

Convergence of higher order finite volume schemes can also be proven within the given 
framework as long as they are L°°-stable and allow for an estimate on the entropy residual in 
the sense of Theorem 2.14, for details see Kroner, Noelle and Rokyta, 1995; Chainais-Hillairet, 
2000. 


2.2.3. Error estimates and convergence rates. There are two primary approaches taken to 
obtain error estimates for approximations of scalar nonlinear conservation laws. One approach 
is based on the ideas of Oleinik and is applicable only in one space dimension (see Oleinik, 
1963; Tadmor, 1991). The second approach which is widely used in the numerical analysis of 
conservation laws is based on the doubling of variables technique of Kruzkov (see Kruzkov, 
1970; Kuznetsov, 1976). In essence, this technique enables one to estimate the error between 
the exact and approximate solution of a conservation law in terms of the entropy residual 
introduced in (44). Thus, an a posteriori error estimate is obtained. Using a 
priori estimates of the approximate solution (see Section 2.2.1, and Theorems 2.9, 2.10), 
a convergence rate or an a priori error estimate is then obtained. The next theorem gives 
a fundamental error estimate for conservation laws independent of the particular finite 
volume scheme (see Eymard, Galluoet and Herbin, 2000; Chainais-Hillairet, 1999; Kroner and 
Ohlberger, 2000). 


Theorem 2.16 (Fundamental error estimate) Let uq £ -£?F(M d ) and let u be an entropy 
weak solution of (la- lb). Furthermore, let v £ L°°(R d x R+ ) be a solution of the following 
entropy inequalities with residual term R: 


[ [ VK(v)d t <f>A-F Vx (v) ■ V^+ [ 'n K (uo)<t>(',0)>~R{K,v ) 

jRd JR+ JR* 


0 ) 


(45) 


for all k £ R and (j) £ Cq (M d x R+ , R + ). Suppose that there exist measures p v £ Al(R d x R + ) 
and v v £ Al(R d ) such that R(k,v,^>) can be estimated independently of k by 


R(k, v , 4>) < {\dt4 > | + | V0I, p v ) + <|0(*, 0)1, v v ). (46) 

Let K CC R d x R + , oj = Lip(f) , and choose T,R and xo such that T E]0,~[ and K lies 
within its cone of dependence Do, i.e. K C Do where Ds is given as 


D s := [J £fl- w t+<5(zo) x {t}. 

o <t<T 


(47) 
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Then , there exists a 5 > 0 and positive constants C \ , C 2 such that u, v satisfy the following 
error estimate 

l|u - v\\ L i( K) < T(i/ v (Br+s{x 0 )) + CiHv{D s ) + C^V^APs)). (48) 

This estimate can be used either as an a posteriori control of the error, as the right-hand side 
of the estimate (48) only depends on v, or it can be used as an a priori error bound if one is 
able to estimate further the measures p v and v v using some a priori bounds on v. Finally, note 
that comparable estimates to (48) are obtainable in an L°°(0, T; L 1 (R d ))-norm (see Cockburn 
and Gau, 1995; Bouchut and Perthame, 1998). 

2.2.4 * A posteriori error estimate. 

Theorem 2.17 (A posteriori error estimate) Assume the conditions and notations as in 
Theorem 2.16. Letv — Uh be a numerical approximation to (la-lb) obtained from a monotone 
finite volume scheme that satisfies a local CFL-like condition as given in (39). Then the 
following error estimate holds 

f \u-U h \ < T(\\u 0 - U/ l (-,0)|| L i(BH + h(xo)) + ClT) + C 2 y/v), (49) 

J K 

where 

E K n+1 - u ?\ Atn ^ + 2^ £ At n (At n + - «P I . 

n€loj<=M(t n ) nelo (j,l)eE(t n ) 

Q, u H 2 jjliW) - 

3 5 ' u — V 

with the index sets 7 0) M(t), E(t) given by 

Jo =: {n I 0<t n <min{^-^,r} , 

LJ 

M(t) H {j | there exists x € Tj such that (x,t) € Dr+s} , 

E(t) = {(jj) | there exists x 6 Tj U Ti such that (x,t) £ D R +s} • 

Furthermore , the constants Ci,C 2 only depend on T } \\uq\\bv ond ||uo|jx°° (for details see 
Kroner and Ohlberger, 2000). 

Note that this a posteriori error estimate is local, since the error on a compact set K is 
estimated by discrete quantities that are supported in the cone of dependence Dr+$. 

2.2.5. A priori error estimate. Using the weak BV estimate (Theorem 2.9) and the Lipschitz 
estimate in time (Theorem 2.10), the right-hand side of the a posteriori error estimate (49) can 
be further estimated. This yields an a priori error estimate as stated in the following theorem 
(for details see Eymard, Galluoet and Herbin, 2000; Chainais-Hillairet, 1999). 

Theorem 2.18 (A priori error estimate) Assume the conditions and notations as in 
Theorem 2.16 and let v = Uh be the approximation to (la), (lb) given by a monotone finite 
volume scheme that satisfies a local CFL-like condition as given in (39). Then there exists a 
constant C > 0 such that 

[ | u-Uhl < Ch 1/4 . 

Jk 

Moreover , in the one- dimensional case , the optimal convergence rate of h 1 / 2 is obtained. 

Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. 

© 2004 John Wiley & Sons, Ltd. 



18 


ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS 


2.2.6 . Convergence proofs via the streamline diffusion discontinuous Galerkin finite element 
method. It is straightforward to show that the fully-discrete finite volume scheme (22) can 
be viewed as a specific case of the more general streamline diffusion discontinuous Galerkin 
(SD-DG) finite element method which utilizes the mesh dependent broken space V* defined 
as 

V* = {v\v\ tGP p (T), VTeT} (50) 

with V P {T) the space of polynomials of degree < p in the control volume T. By generalizing the 
notion of gradient and flux to include the time coordinate as well, the discontinuous Galerkin 
finite element method for a space-time tessellation T n spanning the time slab increment 
[£ n ,t n+1 ] is given compactly by the following variational statement. 

SD-DG(p) finite element method. Find Uh € VJf such that Vvh E Vjf and n = 0, 1, . . . 

52 ( f (Vh + 6(u h )f(u h ) -Vv h ) V * f(u h )dx + I e(uh)Vu h -Vv h dx 

t7tA Jt Jt 

+ v h -(g(u h -,u h+ )-f(u h -)-i')ds)=0 (51) 

J dT ) 

where in the integration over dT it is understood for the portion x E dT D dV that Uh- r Vh- 
denotes the trace restriction of Uh(T) and Vh(T) onto dT and Uh+ denotes the trace restriction 
of Uh(T f ) onto dT Given this space-time formulation, convergence results for a scalar nonlinear 
conservation law in multi-dimensions and unstructured meshes are given in Jaffre, Johnson 
and Szepessy, 1995 for specific choices of the stabilization functions 6(uh) : i ^ 1 + and 
€(u/*) : M 14 M + together with a monotone numerical flux function g{uh- : u/i+). Using their 
stabilization functions together with a monotone flux function, the following convergence result 
is obtained: 

Theorem 2.19 (SD-DG(p) convergence) Suppose that components of f f (u) E C d (M) are 
bounded and that uq E L 2 (M d ) has compact support. Then the solution Uh of the SD-DG(p) 
method converges strongly in Lp OC (R d x Rff), 1 < p < 2, to the unique solution u of the scalar 
nonlinear conservation law system (la-lb) as H = max(||/i||x 00 (Rd), At) tends to zero. 

The proof of convergence to a unique entropy solution on general meshes for p > 0 is 
based on an extension by Szepessy, 1989 of a uniqueness result by DiPerna, 1985 by providing 
convergence for a sequence of approximations satisfying: 

• a uniform Loo bound in time and L 2 in space, 

• entropy consistency and inequality for all Kruzkov entropies, 

• consistency with initial data. 

By choosing SD-DG(O), the dependence on the as yet unspecified stabilization functions 
S(uh ) and e(uh) vanishes identically and the fully-discrete scheme (22) with monotone flux 
function is exactly reproduced, thus yielding a convergence proof for general scalar conservation 
laws for the finite volume method as well. 
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3. Higher order accurate FV generalizations 

Although an OQi 1 / 2 ) Li-norm error bound for the monotone and E-flux schemes of Sect. 2 is 
known to be sharp (Peterson, 1991), an 0(h) solution error is routinely observed in numerical 
experiments with convex flux functions. Even so, first order accurate schemes are generally 
considered too inaccurate for most quantitative calculations unless the mesh spacing is made 
excessively small thus rendering the schemes inefficient. Godunov, 1959 has shown that all 
linear schemes that preserve solution monotonicity are at most first order accurate. The low 
order accuracy of these monotonicity preserving linear schemes has motivated the development 
of higher order accurate schemes with the important distinction that these new schemes utilize 
essential nonlinearity so that monotone resolution of discontinuities and high order accuracy 
away from discontinuities are simultaneously attained. 

3.1. Higher order accurate FV schemes in 1-D 

A significant step forward in the generalization of Godunov’s finite volume method to higher 
order accuracy is due to van Leer, 1979. In the context of Lagrangian hydrodynamics with 
Eulerian remapping, van Leer generalized Godunov’s method by employing linear solution 
reconstruction in each cell (see Fig. 3b). Let N denote the number of control volume cells in 



a. Cell averaging of quartic data b. Linear reconstruction 


c. Quadratic reconstruction 


Figure 3. Piecewise polynomial approximation used in the finite volume method: (a) cell 
averaging of analytic data, (b) piecewise linear reconstruction from cell averages and (c) 
piecewise quadratic reconstruction from cell averages. 

space so that the j- th cell extends over the interval Tj = [xj-i/ 2 i x j+i/ 2 ] with length Ax, such 
that Ui <j<nTj = [0, 1] with T{ HTj = 0, z ^ j. In a purely Eulerian setting, the higher order 
accurate schemes of van Leer are of the form 

where g(u,v) is a numerical flux function utilizing states u” ±1 ^ 2 and ut jtl y 2 obtained from 
evaluation of the linear solution reconstructions from the left and right cells surrounding the 
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interfaces Xj± i/ 2 , By altering the slope of the linear reconstruction in cells, non-oscillatory 
resolution of discontinuities can be obtained. Note that although obtaining the exact solution 
of the scalar nonlinear conservation law with linear initial data is a formidable task, the solution 
at each cell interface location for small enough time is the same as the solution of the Riemann 
problem with piecewise constant data equal to the linear solution approximation evaluated 
at the same interface location. Consequently, the numerical flux functions used in Sect. 2 can 
be once again used in the generalized schemes of van Leer. This single observation greatly 
simplifies the construction of higher order accurate generalizations of Godunov’s method. The 
ideas of van Leer have been extended to quadratic approximations in each cell (see Fig. 3c) by 
Colella and Woodward, 1984. Although these generalizations of Godunov’s method and further 
generalizations given later can be interpreted in 1-D as finite difference discretizations, concepts 
originally developed in 1-D such as solution monotonicity, positive coefficient discretization and 
discrete maximum principle analysis are often used in the design of finite volume methods in 
multiple space dimensions and on unstructured meshes where finite difference discretization is 
problematic. 

3.1.1. TVD schemes. In considering the scalar nonlinear conservation law (la-lb), Lax, 1973 
made the following basic observation: 

“the total increasing and decreasing variations of a differentiable solution between 
any pair of characteristics are conserved”. 

Furthermore, in the presence of shock wave discontinuities, information is lost and the total 
variation decreases. For the 1-D nonlinear conservation law with compactly supported (or 
periodic) solution data u(x,t), integrating along the constant time spatial coordinate at times 
ti and t 2 yields 

/ OO rO O 

\du(x,t 2 )\ < / \du(x,ti)\, t 2 >h . (52) 

-oo J — oo 

This motivated Harten, 1983 to consider the discrete total variation 

TV(Wfc) ^ \^j+l/2 U h\ > -Aj-f- 1/2^6 = u j+l ~ U j 

j 

and the discrete total variation non-increasing (TVNI) bound counterpart to (52) 

TV« +1 ) < TV«) (53) 

in the design of numerical discretizations for nonlinear conservation laws. A number of simple 
results relating TVNI schemes and monotone schemes follow from simple analysis. 

Theorem 3.1 (TVNI and monotone scheme properties, Harten, 1983) (i) Monotone 
schemes are TVNI. (ii) TVNI schemes are monotonicity preserving, i.e. the number of solution 
extrema is preserved in time. 

Property (i) follows from the L x contraction property of monotone schemes. Property (ii) is 
readily shown using a proof by contradiction by assuming a TVNI scheme with monotone 
initial data that produces new solution data at a later time with interior solution extrema 
present. Using the notion of discrete total variation, Harten, 1983 then constructed sufficient 
algebraic conditions for achieving the TVNI inequality (53). 
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Theorem 3.2 (Harten’s explicit TVD criteria) The fully discrete explicit 1-D scheme 
u? +1 = + At (C i+1/2 «) A j+1/2 u n h + D j+ 1 /2 «) Aj-i/ 2 u%) , j = (54) 

is total variation non-increasing if for each j 

c j+ 1/2 > o , (55a) 

^i+i/2 < 0 , (55b) 

1 - At (Cj-i /2 - D j+ 1/2 ) > 0 . (55c) 

Note that although the inequality constraints (55a-55c) in Theorem 3.2 insure that the total 
variation is non-increasing, these conditions are often referred to as total variation diminishing 
(TVD) conditions. Also note that inequality (55c) implies a CFL-like time step restriction that 
may be more restrictive than the time step required for stability of the numerical method. The 
TVD conditions are easily generalized to wider support stencils written in incremental form, 
see for example Jameson and Lax, 1986 and their corrected result in Jameson and Lax, 1987. 

While this simple Euler explicit integration scheme may seem too crude for applications 
requiring true high order space-time accuracy, special attention and analysis is given to this 
fully-discrete form because it serves as a fundamental building block for an important class 
of high order accurate Runge-Kutta time integration techniques discussed in Sect. 4.1 that, 
by construction, inherit TVD (and later maximum principle) properties of the fully-discrete 
scheme. 

Theorem 3.3 (Generalized explicit TVD criteria) The fully discrete explicit 1-D 


scheme 

k-i 

u^ +1 =u^ + AtJ2C^ 1/2 (u n h )A j+l+1/2 ul, j = 1, . . . ,N (56) 

l=-k 

with stencil width parameter k is total variation non-increasing if for each j 

C^-jl > o , (57a) 

cj+ijt < 0, (57b) 

C j+ 1 % - Cj- 1/2 > 0 , —k + l<l<k — 1, 1^0, (57c) 

1 - A t (cj° } 1/2 - C ( } ~1) 2 ) > 0 . (57d) 


The extension to implicit methods follows immediately upon rewriting the implicit scheme 
in terms of the solution spatial increments A j + j +1 / 2 wji and imposing sufficient algebraic 
conditions such that the implicit matrix acting on spatial increments has a nonnegative inverse. 

Theorem 3.4 (Generalized implicit TVD criteria) The fully discrete implicit 1-D 
scheme 

< +1 -A t E C ]+i /2 K +1 )A j+ / +1/ 2< +1 = u", j = l,...,N (58) 

l= — k 
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with stencil width parameter k is total variation non-increasing if for each j 


i) 

1/2 

J+l/2 


S+l/2 S-l/2 


> 

0 

< 

0 

> 

0 


-* + !</<*-!, / 7 ^ 0 . 


(59a) 

(59b) 

(59c) 


Theorems 3.3 and 3.4 provide sufficient conditions for non-increasing total variation of explicit 
(56) or implicit (58) numerical schemes written in incremental form. These incremental 
forms do not imply discrete conservation unless additional constraints are imposed on the 
discretizations. A sufficient condition for discrete conservation of the discretizations (56) and 
(58) is that these discretizations can be written in a finite volume flux balance form 

fc-i 

9j+l/2 - dj-l/2 = ^2 Cfl l / 2 { u h)kj+l+l/2'U'h 
i=-k 

where gj±i /2 are the usual numerical flux functions. Section 3.1.2 provides an example of how 
the discrete TVD conditions and discrete conservation can be simultaneously achieved. A more 
comprehensive overview of finite volume numerical methods based on TVD constructions can 
be found the books by Godlewski and Raviart, 1991 and LeVeque, 2002. 


3.1.2. MUSCL schemes. A general family of TVD discretizations with 5-point stencil is the 
Monotone Upstream-centered Scheme for Conservation Laws (MUSCL) discretization of van 
Leer, 1979; van Leer, 1985. MUSCL schemes utilize a /c-parameter family of interpolation 
formulas with limiter function \&( R ) : E M 

u jj r ij2 = u j "f ^ ^ (l/i?i)Aj + i/2?i/ l H — $ (Rj)Aj_i/2Uh 

u t+i /2 = u i~ — j+i/ 2 Uh - ^-j^V{l/R j+ i)A j+3/2 u h (60) 


where Rj is a ratio of successive solution increments 

o _ A i+l/2Uft 
— 

Aj- l/2«/» 


(61) 


The technique of incorporating limiter functions to obtain non-oscillatory resolution of 
discontinuities and steep gradients dates back to Boris and Book, 1973. For convenience, 
the interpolation formulas (60) have been written for a uniformly spaced mesh although the 
extension to irregular mesh spacing is straightforward. The unlimited form of this interpolation 
is obtained by setting 4>( R ) = 1. In this unlimited case, the truncation error for the 
conservation law divergence in (la) is given by 


Truncation Error = 



(Ax)2 J^ /(u) 


This equation reveals that for k = 1/3, the 1-D MUSCL formula yields an overall spatial 
discretization with 0(Aa; 3 ) truncation error. Using the MUSCL interpolation formulas given 
in (60), sufficient conditions for the discrete TVD property are easily obtained. 
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Theorem 3.5 (MUSCL TVD scheme) The fully discrete 1-D scheme 

•T = «J - <*-./! -»?-.*). i = i » 

with monotone Lipschitz continuous numerical flux function 

9j+ 1/2 = #( U j+l/2’ U J+l/2) 

utilizing the K-parameter family of MUSCL interpolation formulas (60) and (61) is total 
variation non-increasing if there exists a #( R ) such that Vi? € K 

l + a)L^ (62a) 

x — fb 1 — K 

and 


0 < 


1 - K 

m) 


R ~ 


< 2 + a 


(62b) 


with a G [—2, 2 (1 — k)/(1 + k)] under the time step restriction 

At 2 — (2 + a)/c ] dg 


1- 


Axj 


1 — K 


\ du 


> 0 


where 


dg max 
du j 


sup 

S€ t“7-I/2 '“7+1/21 

S «l”7-l/2'“7+l/2l 




Qq \ 

“ g~+( u j-i/2’ u )j ■ 


For accuracy considerations away from extrema, it is desirable that the unlimited form 
of the discretization is obtained. Consequently, the constraint $(1) = 1 is also imposed 
upon the limiter function. This constraint together with the algebraic conditions (62a-b) are 
readily achieved using the well known MinMod limiter, with compression parameter /? 

determined from the TVD analysis 

\£ MM (i?) = max(0, min(i?, /?)) , ft € [1, (3 -«)/(! — /c)] . 


Table I. Members of the MUSCL TVD family of schemes. 


K 

Unlimited Scheme 

Anax 

Truncation Error 

1/3 

-1 

0 

1/2 

Third-Order 
Fully Upwind 
Fromm’s 

Low Truncation Error 

4 
2 
3 

5 

0 

n(A x) 2 £sf(u) 
~^(Az) 2 ;|£r/(u) 


Table I summarizes the MUSCL scheme and maximum compression parameter for a number 
of familiar discretizations. Another limiter due to van Leer that meets the technical conditions 
of Theorem 3.5 and also satisfies $(1) = 1 is given by 

i?+|i?| 


$ VL (i?) = 


1 + |i?| 


Encyclopedia of Computational Mechanics. Edited by Erwin Stein, Rene de Borst and Thomas J.R. Hughes. 
© 2004 John Wiley &c Sons, Ltd. 




24 


ENCYCLOPEDIA OF COMPUTATIONAL MECHANICS 


This limiter exhibits differentiability away from R = 0 which improves the iterative convergence 
to steady state for many algorithms. Numerous other limiter functions are considered and 
analyzed in Sweby, 1984. 

Unfortunately, TVD schemes locally degenerate to piecewise constant approximations at 
smooth extrema which locally degrades the accuracy. This is an unavoidable consequence of 
the strict TVD condition. 

Theorem 3.6 (TVD critical point accuracy, Osher, 1984) The TVD discretizations 
(54), (56) and (58) all reduce to at most first order accuracy at non-sonic critical points , 
i.e. points u* at which f'(u *) ^ 0 and u* x = 0. 

3.1.3. ENO/WENO schemes. To circumvent the degradation in accuracy of TVD schemes 
at critical points, weaker constraints on the solution total variation were devised. To this 
end, Harten proposed the following abstract framework for generalized Godunov schemes in 
operator composition form (see Harten et al., 1986; Harten et al., 1987; Harten, 1989) 

u n h +1 =A-E(r)-R° p (-,u n h ) . (63) 

In this equation, £ VJP denotes the global space of piecewise constant cell averages as defined 
in (20), Rp(x) is a reconstruction operator which produces a cell- wise discontinuous p-th order 
polynomial approximation from the given solution cell averages, E(r) is the evolution operator 
for the PDE (including boundary conditions), and A is the cell averaging operator. Since A 
is a nonnegative operator and E(r) represents exact evolution in the small, the control of 
solution oscillations and Gibbs-like phenomena is linked directly to oscillation properties of 
the reconstruction operator, Rp(x). One has formally in one space dimension 

TV« +1 ) = TV (A • E(t) ■ R° p (-, «£)) < TV(R° p (x; u n h )) 

so that the total variation depends entirely upon properties of the reconstruction operator 
Rp(x'jU%). The requirements of high order accuracy for smooth solutions and discrete 
conservation give rise to the following additional design criterion for the reconstruction operator 

• Rp(x; Uh) = u(x) + e(x) Ax p+1 + 0( Ax p+2 ) whenever u is smooth (64a) 

• A\TjRp(x ; uh) = Uh | Tj = Uj, j = 1, . . . , N to insure discrete conservation (64b) 

• TV(R(x;Uh)) < TV(u£) + 0( Ax p+1 ) an essentially non-oscillatory reconstruction. (64c) 

Note that e(x) may not be Lipschitz continuous at certain points so that the cumulative error 
in the scheme is 0(Ax p ) in a maximum norm but remains 0( Ax p+1 ) in an Li-norm. To 
achieve the requirements of (64a-64c), Harten and coworkers considered breaking the task into 
two parts 

• Polynomial reconstruction from a given stencil of cell averages 

• Construction of a “smoothest” polynomial approximation by a solution adaptive stencil 
selection algorithm. 

In the next section, a commonly used reconstruction technique from cell averages is considered. 
This is then followed by a description of the solution adaptive stencil algorithm proposed by 
Harten et al., 1986. 
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3. 1-4- Reconstruction via primitive function. Given cell averages Uj of a piecewise smooth 
function u(x ), one can inexpensively evaluate pointwise values of the primitive function U(x) 

U(x)= j u(0d4 

3 XQ 

by exploiting the relationship 

3 

AxjUj = U(x j+1/2 ) ■ 

3=30 

Let H p (x;u) denote a p-th order piecewise polynomial interpolant of a function u. Since 

u(x) = , 

an interpolant of the primitive function given pointwise samples U(xj+ 1 / 2 ) yields a 
reconstruction operator 

s R 0 p (x;u h ) = ±H p+1 (x;U) . 

As a polynomial approximation problem, whenever U (x) is smooth one obtains 

f^H p (x; U) = f^U(x) + 0(Ax p+1 ~ k ) ,0<k<p 

and consequently 

^R P (x-,u h ) = ^ u(x ) + 0(Ax p +'-‘) . 

By virtue of the use of the primitive function U (x), it follows that 

A\ Tj R Q p (x]u h ) = Uj 

and from the polynomial interpolation problem for smooth data 

R° p (x-,u h ) = u(x) + 0(Ax p+1 ) 

as desired. 

3.1.5. ENO reconstruction. The reconstruction technique outlined in Section 3.1.4 does not 
satisfy the oscillation requirement given in (64c). This motivated Harten and coworkers 
to consider a new algorithm for essentially non-oscillatory (ENO) piecewise polynomial 
interpolation. When combined with the reconstruction technique of Section 3.1.4, the resulting 
reconstruction then satisfies (64a-c). Specifically, a new interpolant H p (x ; u) is constructed so 
that when applied to piecewise smooth data v(x) gives high order accuracy 

v) = f^v(x) + 0(Ax p+1 ~ k ) , 0 <k<p 
but avoids having Gibbs oscillations at discontinuities in the sense 

TV(H p {x-,v)) < TV{v) + 0( Ax p+1 ) . 
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The strategy pursued by Harten and coworkers was to construct such an ENO polynomial 
H p (x]w) using the following steps. Define 

H™°{x;w) = Ppj+i/ 2 (x;w) for x, < x < x j+ i , j = l,...,N 

where P^^/ 2 the pth degree polynomial which interpolates w{x) at the p + 1 successive 
points { Xi }, i p (j ) < i < i p (j) + p that include xj and Xj+ i, i.e. 

Ppj +1 / 2 (*iJ w) = w(xi) , i p (j) < i < ip(j) + p , 1 -p< ip(j) - j < o . (65) 

Equation (65) describes p possible polynomials depending on the choice of i p (j) for an interval 
(xj,Xj+ 1 ). The ENO strategy selects the value i p (j) for each interval that produces the 
“smoothest” polynomial interpolant for a given input data. More precisely, information about 
smoothness of w(x ) is extracted from a table of divided differences of w(x) defined recursively 
for i = 1, . . . , N by 

w(Xi) 

Xi-\-l Xi 

W [xj-j-i , » . • , w[Xi , . . . , Xi^.^—1 j 

Xi+k 

The stencil producing the smoothest interpolant is then chosen hierarchically by setting 

h(j) =j 


U}[Xi] = 
w[xi,x i+1 ] = 

w[Xi,..., x i+k ] = 


and for 1 < k < p — 1 

,• (;\ _ - 1 if ffifi'l 

k+l[]) ~ \ i k (j) otherwise. W 

Harten et al., 1986 demonstrate the following properties of this ENO interpolation strategy 

• The accuracy condition 

P Id+l/ 2 ( X ) = W ( X ) + 0(Ax P+1 ) , X <= (Xj,Xj+l) . 

• Pp NO (x) is monotone in any cell interval containing a discontinuity. 

• There exists a function z(x ) nearby P^ NO (x) in the interval (xj,Xj+ 1 ) in the sense 

z(x) = Ppj+ 1/2 (x) + 0(Ax p+1 ) , x 6 (Xj,x j+ i) 
that is total variation bounded, i.e. the nearby function z[x) satisfies 

TV(z) < TV{w) . 

3.1.6 . WENO reconstruction. The solution adaptive nature of the ENO stencil selection 
algorithm (66) yields non-differentiable fluxes that impede convergence to steady state. In 
addition, the stencil selection algorithm chooses only one of p possible stencils and other 
slightly less smooth stencils may give similar accuracy. When w(x) is smooth, using a linear 
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combination of all p stencils with optimized weights yields a more accurate 0(Ax 2p ~ l ) 
interpolant. More specifically, let P p k j + i/ 2 denote the unique pdlynomial interpolating p + 1 
points with stencil {xj+i- p +k,Xj+i+k} then 

p - i v - i 

P pJ+i/2 (w(x)) !/»(*) + ^Az 2 *’- 1 ) , 53a;* = 1 . 

k—Q k = 0 

For example, optimized weights for p = 1,2,3 yielding 0(Aa; 2p ” 1 ) accuracy are readily 
computed 


P= 1 : 

Wo = 

1, 





2 

1 


P—2: 

Wo = 

3 , *1 = 

3’ 




3 

3 

1 

P = 3: 

w 0 = 

To’ “ 1 = 

= 5’ 

s 

II 

SI 


In the WENO schemes of Jiang and Shu, 1996; Shu, 1999, approximate weights, w, are devised 
such that for smooth solutions 

w* = w* 4- 0( Ax p ~ l ) 

so that the 0{Ax 2p ~' 1 ) accuracy is still retained using these approximations 


p P ,j+ 1/2M*)) = 53 7+1/2 ( x ) + ^(Ax 2p *) , 53^* = 1 • 

k -0 &=0 

The approximate weights are constructed using the ad hoc formulas 

w* ~ 


<*k = 


(e + A) 


2 ’ 


W* = 


ELo a k 


where e is an approximation to the. square root of the machine precision and 0k is a smoothness 
indicator 


0k 


r*;+ 1/2 
/=1 '* X 3 - 1/2 


Axf" 1 




dx . 


For a sequence of smooth solutions with decreasing smoothness indicator /?*, these formulas 
approach the optimized weights, ► cjk- These formulas also yield vanishing weights 

u)k — ► 0 for stencils with large values of the smoothness indicator such as those encountered 
at discontinuities. In this way, the WENO construction retains some of the attributes of the 
original ENO formulation but with increased accuracy in smooth solution regions and improved 
differentiability often yielding superior robustness for steady state calculations. 


3.2. Higher order accurate FV schemes in multi- dimensions. 

Although the one-dimensional TVD operators may be readily applied in multi-dimensions on 
a dimension-by-dimension basis, a result of Goodman and LeVeque, 1985 shows that TVD 
schemes in two or more space dimensions are only first order accurate. 
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Theorem 3.7 (Accuracy of TVD schemes in multi-dimensions) Any two-dimensional 
finite volume scheme of the form 

~ “ 17 T. . (9i+i/2,j - 9i- 1/2, j) “ f-i /2 “ ^i t j -1/2) j 1 < z ^ 1 < J < N 


|r|ij +1/2,<? |T|i j J - 

with Lipschitz continuous numerical fluxes for integers p , 9 , r, 5 

9i+lj2 f j = 9{ u i-p,j-g? ' - * i u i+r t j+s)) 

hiJ+1/2 ~ h{ u i—p,j—qi • * * J ^t+rj+s)? 

that is total variation non-increasing in the sense 

TV{ul +1 ) < TV{ul) 

where 

TV(u) =z [Ay^ + iy 2 ,j Wi+lj — u i,j\ "b \ u iJ+l ~~ u 


1 < i < M, 1 < j < N 


is at most first-order accurate. 

Motivated by the negative results of Goodman and LeVeque, weaker conditions yielding 
solution monotonicity preservation have been developed from discrete maximum principle 
analysis. These alternative constructions have the positive attribute that they extend to 
unstructured meshes as well. 

3.2.1. Positive coefficient schemes on structured meshes. Theorem 2.6 considers schemes of 
the form 

u? + 1 =u? + -£- E <?;*«)« -«") , VTj € T 

1 jl ve^sar,- 

and provides a local space-time discrete maximum principle 

s u >” +1 4 v.&W’O 

VTj E T under a CFL-like condition on the time step parameter if all coefficients Cjk are 
nonnegative. Schemes of this type are often called positive coefficient schemes or more simply 
positive schemes. To circumvent the negative result of Theorem 3.7, Spekreijse, 1987 developed 
a family of high order accurate positive coefficient schemes on two-dimensional structured 
M x N meshes. For purposes of positivity analysis, these schemes are written in incremental 
form on a M x N logically rectangular 2-D mesh 

^ ( Afoj « +li , - ulj) + B ? J+1 «, +1 - 
+ Cl l lti K_ hj - u?j) + D?j_ 1 - «?j) ) , 1 < < < M, 1 < i < N 


where the coefficients are nonlinear functions of the solution 


A n 

l,j 

= A(... 

»«?-i 


Tjn 

&i,j + 1 

= £(.. 

5 ^i,j— 1 5 U i,j > > 


fin 

1J 

= C(... 


...) 

D n 

U ij - 1 

= D(... 

j — 1 ’ » ^zj+1 » 

...) 


Once written in incremental form, the following theorem follows from standard positive 
coefficient maximum principle analysis. 
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Theorem 3.8 (Positive coefficient schemes in multi-dimensions) The discretization (67) 
is a positive coefficient scheme if for each 1 < i < M, 1 < j < N and time slab increment 
[t n ,f n+1 ] 

> 0, B? J+1 > 0, CILu > 0, > 0, (68) 

and 

1 - At {A? +1J + Bl J+1 + C^j + 2?^_j) > 0 (69) 

with discrete space-time maximum principle 




and discrete maximum principle at steady state 


where u * denotes the numerical steady state . 

Using a procedure similar to that used in the development of MUSCL TVD schemes in 1-D, 
Spekreijse, 1987 developed a family of monotonicity preserving MUSCL approximations in 
multi-dimensoins from the positivity conditions of Theorem 3.8. 


Theorem 3.9 (MUSCL positive coefficient scheme) The fully discrete 2-D finite vol- 
ume scheme 


u i,j \T\ij | T |. j (Kj+i /2 1 / 2 ) » 

utilizing monotone Lipschitz continuous numerical flux functions 

9i+l/2,j — 

h i,j+ 1/2 “ h{ U i y j+i/2' U Xj+l/2) 

and MUSCL extrapolation formulas 


1 <i<M, 1 <j<N 


U i+l/2,j 

U t+l/ 2 ,j 

U i,j+l/2 

u tj+ 1/2 


u hj 2 ( u hj ~ 

u hj 'h 

u i,j i U i,j ~ u i,j— 1 ) 

u *,j 2 V ^»,i) (^J “ u *J+l) 


where 


c> — 
riij — 


J i,j 


'U'i,j + 1 j 


satisfies the local maximum principle properties of Lemma 3.8 and is second order accurate if 
the limiter = 'T ( /?) ftas f/ie properties that there exist constants 3 € (0, oo), a € [—2,0] suc/i 
that ViZ € M 


a < *(/?) < 0 , 


~ 0 < 


m) 

R 


<2 + a 


(70) 
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with the constraint \t(l) = 1 and the smoothness condition \P( R ) 6 C 2 near R = 1 together 
with a time step restriction for stability 


where 

max 

i,j 

max 
hj 


££ 

dh 

du 


i — (i +<&) 


At 


d£ 

<9u 


+ 




dh 

du 


n,max 

Lj 


>0 




sup 

h u + i 


sup 


2€[u". 


,j-l/2’“i,i + l/2 J 


d~=^ U t+ 1/2, j) au +K-l/2,j> W )) ^ 0 


dh + . dh . xV 

^r( U - U i, J+ l/2) 5 u+ Kj-i/ 2> u )J ^ 




Many limiter functions satisfy the technical conditions (70) of Theorem 3.9. Some examples 
include 


• the van Leer limiter 


• the van Albada limiter 


«f VL (i?) 


R + \R\ 
1 + 1*1 


tf VA (H) 


R + R 2 
l + R 2 ' 


In addition, Koren, 1988 has constructed the limiter 


# K (i?) 


R + 2R 2 
2 - R + 2R? 


which also satisfies the technical conditions (70) and corresponds for smooth solutions in 1-D 
to the most accurate n = 1/3 MUSCL scheme of van Leer. 


3.2.2. FV schemes on unstructured meshes utilizing linear reconstruction. Higher order finite 
volume extensions of Godunov discretization to unstructured meshes are of the general form 

^ = -pr E • VT ^ T < 71 > 

with the numerical flux gjk{u,v) given by the quadrature rule 

Q 

9jk( u jki u ^k) = UJq SiVjki^q)] u jk( x q)> u ^k( x <i)) > 

9 =i 

where € R and € ejk represent quadrature weights and locations, q = 1, . . . , Q. Given the 
global space of piecewise constant cell averages, u/i 6 Vjf, the extrapolated states 1 /^( 2 ) and 
u+ k (x) are evaluated using ap-th order polynomial reconstruction operator, R® : *-* VJf, 

u~ k (x) = lim i?p (a; - e (x) ;u h ) 
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Figure 4. Polygonal control volume cell Tj and 
perimeter quadrature points (solid circles). 

uf k (x) = lim-Rp(x + e v jk (x)\u h ) 

for x € ejk- In addition, it is assumed that the reconstruction satisfies the property 
]T~l It Rp( x ’ u h)dx = Uj stated previously in (64b). In the general finite volume formulation, 
the control volume shapes need not be convex, see for example Fig. 4. Even so, the solution 
accuracy and maximum stable time step for explicit schemes may depend strongly on the 
shape of individual control volumes. In the special case of linear reconstruction, i?? (x; uh ), the 
impact of control volume shape on stability of the scheme can be quantified more precisely. 
Specifically, the maximum principle analysis presented later for the scheme (71) reveals an 
explicit dependence on the geometrical shape parameter 

T° eom = sup a~ l (9) (73) 

O<0<2 7T 

where 0 < a(6) < 1 represents the smallest fractional perpendicular distance from the 
gravity center to one of two minimally separated parallel hyperplanes with orientation 9 and 
hyperplane location such that all quadrature points in the control volume lie between or on the 
hyperplanes as shown in Fig. 5. Table II lists T 9eom values for various control volume shapes in 
R 1 , K 2 , R 3 , and R rf . As might be expected, those geometries that have exact quadrature point 
symmetry with respect to the control volume gravity center have geometric shape parameters 
pgeom e q ua i 2 regardless of the number of space dimensions involved. 

Lemma 3.10 (Finite volume interval bounds on unstructured meshes, R,5( x ; u *0) 
The fully discrete finite volume scheme 

U T 1=U ?-W\ £ 9Au- k \u^ n ) , VI} € T (74) 

1 31 Vejfc <EdTj 

with monotone Lipschitz continuous numerical flux function , nonnegative quadrature weights , 
and linear reconstructions 

Uj k {x) 

u jk( x ) 
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£ gravity center 
• quadrature point 


Figure 5. Minimally separated hyperplanes h L (0) and h u (0) and the 
fractional distance ratio a(0) for use in the calculation of r seom . 


Table II. Reconstruction geometry factors for various control 
volume shapes utilizing midpoint quadrature rule. 


control volume shape 

space dimension 

pgeom 

segment 

i 

2 

triangle 

2 

3 

parallelogram 

2 

2 

regular n-gon 

2 

"Wl 

tetrahedron 

3 

4 

parallelepiped 

3 

2 

simplex 

d 

d+1 

hyper-parallelepiped 

d 

2 

polytope 

d 

Eqn. (73) 


with extremal trace values at control volume quadrature points 


U f n = mm r , u%{x q ) , t/f ax = y max, uf k (x q ) , x q € e jk 


i<q<Q 


1 <q<Q 


exhibits the local interpolated interval bound 

ajUf n ' n + (1 - aj)uJ < u] +l < (1 - aj) u" + a 3 U] 
with the time step proportional interpolation parameter aj defined by 

9g 


Vi 


s i^i r ! '°” E 


sup 

min,n ^max,n 


ye jk ea Tj seivf 


du + 


("j k(x q ));u,u) 


that depends on the shape parameter r geom defined in (73). 


(75) 


(76) 
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Given the two-sided bound of Lemma 3.10, a discrete maximum principle is obtained under 
a CFL-like time step restriction if the limits C/j 11 ** and U™ in can be bounded from above and 
below respectively by the neighboring cell averages. This idea is given more precisely in the 
following theorem. 

Theorem 3.11 (Finite volume maximum principle on unstructured meshes, Rj) Let 

u™ m and ii™ ax denote the minimum and maximum value of solution cell averages for a given 
cellTj and corresponding adjacent cell neighbors , i.e. 

uf n = mm (u j ,u k ) and uf™ = max (u jt u k ) . (77) 

Vej k edTj J VejhEoTj 

The fully discrete finite volume scheme 

“< +1 = “”"i Pi £ »*<»«". «j') . W>€ r (78) 

1 jl Ve, k edT, 


with monotone Lipschitz continuous numerical flux function, nonnegative quadrature weights, 
and linear reconstructions 

uj k (x) = limi?£(x -eu jk (x);u h ) , x € e jk , u h €Vg 

u tk( x ) - + 'jk(x);u h ) , x € e jk , u h € (79) 

exhibits the local space-time maximum principle for each Tj E T 

min (uf.ul) < u 7 }* 1 < max (u^uT) 

1 ~ 3 “Ve ifc €Wi- 3 k ' 


as well as the local spatial maximum principle at steady state (u n+1 = u n = u*) 


min ui < u*j < 

Ve^edTj ~ 3 ~ 


max ut 

Ve jk edTf 


if the linear reconstruction satisfies Vej k E dTj and x q E q = 1, . . . , Q 
ma x(iif < uj£ n (z q ) < min (u^ n t u^ n ) 

under the time step restriction 


i _ ^ pgeom 

I Tji ^ 


sup 


v *i* €aT i 
1 <q<Q 


min, n 


x _ f n 


(^*(x„));S,u)| > 0 


m in , n max, n , 


(80) 


with r geom defined in (73). 

Note that a variant of this theorem also holds if the definition of u max and u min are expanded 
to include more control volume neighbors. Two alternative definitions frequently used when 
the control volume shape is a simplex are given by 

uf in = min Uk and idp ax = max u k . (81) 

3 T fc €T 3 T k €T y ' 
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These expanded definitions include adjacent cells whose intersection with Tj in need only 
be a set of measure zero or greater. 

Slope limiters for linear reconstruction. Given a linear reconstruction Ri(x;Uh) that does 
not necessarily satisfy the requirements of Theorem 3.11, it is straightforward to modify 
the reconstruction so that the new modified reconstruction does satisfy the requirements of 
Theorem 3.11. For each control volume Tj eT a modified reconstruction operator R\{x]Uh) 
of the form 

Rl{x-,u h )\ Ti = Uj + aTj (-R? (x; -Uj) 

is assumed for ax, £ [0, 1]. By construction, this modified reconstruction correctly reproduces 
the control volume cell average for all values of a t, , i.e. 

A- [ j?°( X ;u A )dx = Uj . (82) 

Mil Jt 5 

The most restrictive value of ax, for each control volume Tj is then computed based on the 
Theorem 3.11 constraint (80), i.e. 


a£? M = 


mm 

V'jktdTj 
1 <q<Q 


f min(ti^ iax ,ng iax )~Uj 
i ma x(uf in 

RfiXq^Uh^Tj-Uj 


if R%(x q ;u h )\ Tj > min(u^ ax ,u^ ax ) 

if < max(uf in ,uf in ) 

otherwise 


(83) 


where u max and u min are defined in (77). When the resulting modified reconstruction operator 
is used in the extrapolation formulas (79), the discrete maximum principle of Theorem 3.11 is 
attained under a CFL-like time step restriction. By utilizing the inequalities 

max(%,Ujb) < min (u™^,u™ ax ) and min (uj,u k ) > ma x(u“ m ,ti™ in ) 


it is straightforward to construct a simpler but more restrictive limiter function 



mm 

1<<?<Q 


f ma x(uj Uj 

Rl{x q \U h ]\Tj-Uj 
/ min(nj,tife) — Uj 
R°(x q ;u h )\ Tj -Uj 

1 


if #?(x 9 ;u/,)|t, > ma x(uj,u k ) 
if Ri(x q ;uh)\Tj <mm(v,j,u k ) 
otherwise 


(84) 


that yields modified reconstructions satisfying the technical conditions of Theorem 3.11. This 
simplified limiter (84) introduces additional slope reduction when compared to (83). This can 
be detrimental to the overall accuracy of the discretization. The limiter strategy (84) and other 
variants for simplicial control volumes are discussed further in Liu, 1993; Wierse, 1994; Batten, 
Lambert and Causon, 1996. 

In Barth and Jespersen, 1989, a variant of (83) was proposed 


O.T- — min 

3 V.j.fcgW, 

!<4<Q 


r u™* x ~uj 
R'i(x q ;u h )\ Tj -u j 

U™ ln ~Uj 

RU X 9'y U h)\T j ~Uj 


if fJ?(x g ;uh)|Tj > uf ax 

if Ri(x q -,u h )\ Tj < uf n 
otherwise 


(85) 


Although this limiter function does not produce modified reconstructions satisfying the 
requirements of Theorem 3.11, using Lemma 3.10 it can be shown that the Barth and Jespersen 
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limiter yields finite volume schemes (74) possessing a global extremum diminishing property, 

i.e. that the solution maximum is non-increasing and the solution minimum is nondecreasing 
between successive time levels. This limiter function produces the least amount of slope 
reduction when compared to the limiter functions (83) and (84). Note that in practical 
implementation, all three limiters (83), (84) and (85) require some modification to prevent 
near zero division for nearly constant solution data. 

3.2.3. Linear reconstruction operators on simplicial control volumes. Linear reconstruction 
operators on simplicial control volumes that satisfy the cell averaging requirement (64b) often 
exploit the fact that the cell average is also a pointwise value of any valid linear reconstruction 
evaluated at the gravity center of a simplex. This reduces the reconstruction problem to that of 
gradient estimation given pointwise samples at the gravity centers. In this case, it is convenient 
to express the reconstruction in the form 

f^(x; u h ) |t,- = uj + (Vtifc)T, • (x - xj) (86) 

where xj denotes the gravity center for the simplex Tj and (Vu/Jt, is the gradient to be 
determined. Figure 6 depicts a 2-D simplex A123 and three adjacent neighboring simplices. 
Also shown are the corresponding four pointwise solution values {A, B , C, 0} located at gravity 
centers of each simplex. By selecting any three of the four pointwise solution values, a set of 
four possible gradients are uniquely determined, i.e. {V(ABC), V(ABO), V(BCO), V(CAO )}. 
Using the example of Fig. 6, a number of slope limited reconstruction techniques are possible 



1 


Figure 6. Triangle control volume A 123 (shaded) with three adjacent cell neighbors. 

for use in the finite volume scheme (78) that meet the technical conditions of Theorem 3.11. 

1. Choose (Vu/ 1 )t 123 = V(ABC) and limit the resulting reconstruction using (83) or (84). 
This technique is pursued in Barth and Jespersen, 1989 but using the limiter (85) instead. 

2. Limit the reconstructions corresponding to gradients V(ABC'), V(ABO), V(J5CO) and 
V(CAO) using (83) or (84) and choose the limited reconstruction with largest gradient 
magnitude. This technique is a generalization of that described in Batten, Lambert and 
Causon, 1996 wherein limiter (84) is used. 

3. Choose the unlimited reconstruction V(ABC),V{ABO),V(BCO) and V(CAO) with 
largest gradient magnitude that satisfies the maximum principle reconstruction bound 
inequality (80). If all reconstructions fail the bound inequality, the reconstruction 
gradient is set equal to zero, see Liu, 1993. 
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Figure 7. Triangulation of gravity center locations showing a 
typical control volume To associated with the triangulation vertex 
vq with cyclically indexed graph neighbors T*, k = 1, . . . , Nq. 


3.2.4 • Linear reconstruction operators on general control volumes shapes. In the case of linear 
reconstruction on general volume shapes, significant simplification is possible when compared 
to the general p-exact reconstruction formulation given in Sect. 3.2.5. It is again convenient to 
express the reconstruction in the form 

R°i(x; u h )\ Tj = Uj + (Vu h ) T} ■ (X - x®) (87) 

where x ? denotes the gravity center for the control volume Tj and is the gradient to be 

determined. Two common techniques for simplified linear reconstruction include a simplified 
least squares technique and a Green-Gauss integration technique. 

Simplified least squares linear reconstruction. As was exploited in the linear reconstruction 
techniques for simplicial control volumes, linear reconstructions satisfying (64b) on general 
control volume shapes are greatly simplified by exploiting the fact that the cell average value 
is also a pointwise value of all valid linear reconstructions evaluated at the gravity center of 
a general control volume shape. This again reduces the linear reconstruction problem to that 
of gradient estimation given pointwise values. In the simplified least squares reconstruction 
technique, a triangulation (2D) or tetrahedralization (3D) of gravity centers is first constructed 
as shown in Fig. 7. Referring to this figure, for each edge of the simplex mesh incident to the 
vertex v 0 , an edge projected gradient constraint equation is constructed subject to a pre- 
specified nonzero scaling Wk 


Wk (Vti/OTo • (x% - X®) = w k (u k - u 0 ) . 

The number of edges incident to a simplex mesh vertex in R d is greater than or equal to d 
thereby producing the following generally non-square matrix of constraint equations 


w\ Axf wi Ay® 


/ wi(m ~u 0 ) \ 

W No Ax No w No Ay% 0 _ 

(Vu*) To = 

\WN 0 (uNo - «o)/ 


or in abstract form 

[U U\ V« = / . 
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This abstract form can be symbolically solved in a least squares sense using an 
orthogonalization technique yielding the closed form solution 


1 ( *22 (^I ■ f) - 11202 ■ f)\ 

*11*22 - l \ 2 \ln02 ■ f) ~ *12 {01 ■ f) ) 


( 88 ) 


with lij = L{ • Lj. The form of this solution in terms of scalar dot products over incident edges 
suggests that the least squares linear reconstruction can be efficiently computed via an edge 
data structure without the need for storing a non-square matrix. 

Green-Gauss linear reconstruction . This reconstruction technique specific to simplicial meshes 
assumes nodal solution values at vertices of the mesh which uniquely describes a C° linear 
interpolant, Uh> Gradients are then computed from the mean value approximation 

|f*o| (Vu/,)n 0 « f Vu/, dx = t U h dv . (89) 

For linear interpolants, the right-hand side term can be written in the following equivalent 


k+l 



Figure 8. Median dual control volume To demarcated by median segments of triangles 
incident to the vertex vq with cyclically indexed adjacent vertices t/*., k = 1, . . . No. 


form using the configuration depicted in Fig. 8 



• N 0 

Vu h dx = -(uq + u k ) i/ ok 

k - 1 Z 


where v Qk represents any path integrated normal connecting pairwise adjacent simplex gravity 
centers, i.e. 

f X k + 1/2 

vok = / dv . (90) 

A particularly convenient path is one that traces out portions of median segments as shown in 
Fig. 8. These segments demarcate the so-called median dual control volume. By construction, 
the median dual volume |T 0 | is precisely equal to |fi 0 |/3 in 2-D. Consequently, a linear 
reconstruction operator on non-overlapping median dual control volumes is given by 

No 1 

PoKV'U/OTq ~ o(^0 * (91) 

*=1 z 
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The gradient calculation is exact whenever the numerical solution varies linearly over the 
support of the reconstruction. Since mesh vertices are not located at the gravity centers of 
median dual control volumes, the cell averaging property (64b) and the bounds of Theorem 
3.11 are only approximately satisfied using the Green-Gauss technique. 

A number of slope limited linear reconstruction strategies for general control volume shapes 
are possible for use in the finite volume scheme (78) that satisfy the technical conditions of 
Theorem 3.11. Using the example depicted in Fig. 7, let V* + i/ 2 «/i denote the unique linear 
gradient calculated from the cell average set {uo,u*;,Ufc+i}. Three slope limiting strategies 
that are direct counterparts of the simplex control volume case are: 

1. Compute (Vuft)To using the least squares linear reconstruction or any other valid linear 
reconstruction technique and limit the resulting reconstruction using (83) or (84). 

2. Limit the reconstructions corresponding to the gradients V^ + i/ 2 u/ l ,fc = 1 , and 

(Vu h ) To using (83) or (84) and choose the limited reconstruction with largest gradient 
magnitude. 

3. Choose the unlimited reconstruction from Vk+i/ 2 u h, h = and (Vu^)t 0 with 

largest gradient magnitude that satisfies the maximum principle reconstruction bound 
inequality (80). If all reconstructions fail the bound inequality, the reconstruction 
gradient is set equal to zero. 

3.2.5. General p-exact reconstruction operators on unstructured meshes. Abstractly, the 
reconstruction operator serves as a finite-dimensional pseudo inverse of the cell averaging 
operator A whose j-th component Aj computes the cell average of the solution in Tj 

udx . 

The development of a general polynomial reconstruction operator, R that reconstructs p- 
degree polynomials from cell averages on unstructured meshes follows from the application of 
a small number of simple properties. 

1. (Conservation of the mean) Given solution cell averages u h j the reconstruction R p Uh is 
required to have the correct cell average, i.e. 

if v = RpUh then Uh—Av . 

More concisely, 

AR° P = I 

so that R p is a right inverse of the averaging operator A. 

2. (p-exactness) A reconstruction operator is p-exact if R p A reconstructs polynomials 
of degree p or less exactly. Denoting by V v the space of all polynomials of degree p, 

if u eV p and v = An then R° p v ~ u . 

This can be written succinctly as 

R° P A\ Vp = I 

so that Rp is a left inverse of the averaging operator A restricted to the space of 
polynomials of degree at most p. 
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3. (Compact support) The reconstruction in a control volume Tj should only depend of cell 
averages in a relatively small neighborhood surrounding Tj . Recall that a polynomial 
of degree p in contains ( p ^ d ) degrees of freedom. The support set for Tj is required 
to contain at least this number of neighbors. As the support set becomes even larger 
for fixed p, not only does the computational cost increase, but eventually the accuracy 
decreases as less valid data from further away is brought into the calculation. 

Practical implementations of polynomial reconstruction operators fall into two classes: 

• Fixed support stencil reconstructions. These methods choose a fixed support set as 
a preprocessing step. Various limiting strategies are then employed to obtain non- 
oseillatory approximation, see for example Barth and Frederickson, 1990; Delanaye, 1996 
for further details. 

• Adaptive support stencil reconstructions. These ENO-like methods dynamically choose 
reconstruction stencils based on solution smoothness criteria, see for example Harten 
and Chakravarthy, 1991; Vankeirsblick, 1993; Abgrall, 1994; Sonar, 1997; Sonar, 1998 
for further details. 


3.2.6. Positive coefficient schemes on unstructured meshes Several related positive coefficient 
schemes have been proposed on multi-dimensional simplicial meshes based on one-dimensional 
interpolation. The simplest example is the upwind triangle scheme as introduced by Billey 
et al., 1987; Desideri and Dervieux, 1988; Rostand and Stoufflet, 1988 with later improved 
variants given by Jameson, 1993; Cournede, Debiez and Dervieux, 1998. These schemes are 
not Godunov methods in the sense that a single multi-dimensional gradient is not obtained 
in each control volume. The basis for these methods originates from the gradient estimation 
formula (91) generalized to the calculation of flux divergence on a median dual tessellation. 
In deriving this flux divergence formula, the .assumption has been made that flux components 
vary linearly within a simplex yielding the discretization formula 


f div(/) dx= f f -dv = ^ ]- ( f(uj ) + f(u k )) ■ u jk 

JT > JdT i VejkedTj Z 


where Ujk is computed from a median dual tessellation using (90). This discretization is the 
unstructured mesh counterpart of central differencing on a structured mesh. Schemes using 
this discretization of flux divergence lack sufficient stability properties for computing solutions 
of general nonlinear conservation laws. This lack of stability can be overcome by adding 
suitable diffusion terms. One of the simplest modifications is motivated by upwind domain 
of dependence arguments yielding the numerical flux 

Qjk(Uj,Uk) = — (/(Uj) 4* /(Ufc)) * Ujk — — \o\jk AjkU (92) 

with djk a mean value (a.k.a. Murman-Cole) linearization satisfying 


v jk • A jkf = a jk A j k u . 

Away from sonic points where f(u*) ~ 0 for u* £ [v,j,Uj+ 1 ], this numerical flux is formally an 
E-flux satisfying (28). With suitable modifications of aj k near sonic points, it is then possible 
to produce a modified numerical flux that is an E-flux for all data, see Osher, 1984. Theorems 
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2.6, 2.7 and 2.8 show that schemes such as (22) using E-fluxes exhibit local discrete maximum 
principles and L 0 0 stability. 

Unfortunately, schemes based on (92) are too dissipative for most practical calculations. 
The main idea in the upwind triangle scheme is to add anti-diffusion terms to the numerical 
flux function (92) such that the sum total of added diffusion and anti-diffusion terms in the 
numerical flux function vanish entirely whenever the numerical solution varies linearly over the 
support of the flux function. In all remaining situations, the precise amount of anti-diffusion 
is determined from maximum principle analysis. 



Figure 9. Triangle complex used in the upwind triangle schemes showing the linear 
extension of ejk into neighboring triangle for the determination of points xy and x k > . 


Theorem 3.12 (Maximum Principles for the Upwind Triangle Scheme) Let T de- 
note the median dual tessellation of an underlying simplicial mesh. Also let Uj denote the nodal 
solution value at a simplex vertex in one-to-one dual correspondence with the control volume 
Tj € T such that a C° linear solution interpolant is uniquely specified on the simplicial mesh . 
Let gjk(uj',Uj,Uk,Uk>) denote the numerical flux function with limiter function $(•) rl^l 


g jk {u r ,u h u k ,u k ,) = l -{f( Uj ) + f(u k )) ■ v jk - ^ af k (l - A jk u 

+ 5*0 -*(&€?))**•■ 


utilizing the mean value speed ajk satisfying 

Vjk ‘ A jkf = Q>jh 


and variable spacing parameter hjk = |Aj* 2 ;|. The fully discrete finite volume scheme 


uf- 1 = u” 


■pyTT m( U "’> U j> U k,Uk,) 

1 3 } Ve jk edTj 


VTj 6 T 


with linearly interpolated values Uj* and u & as depicted in Fig. 9 exhibits the local space-time 
maximum principle 

min (u^ilZ) < u? +1 < max (u?,u£) 

Ve jk 6dTj 3) ~ 3 ~~ Ve jk edTj 3 

and the local spatial maximum principle at steady state (u n+1 = u n = u*) 


min ut < u* < 
Ve jk edTj * - 3 - 


max 

Ve jh edTj 


Uk 
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if the limiter ^f(R) satisfies Vi? € R 

0 <V(R)/R , 0<$(i?)<2 . 


Some standard limiter functions that satisfy the requirements of Theorem 3.12 include 


• the MinMod limiter with maximum compression parameter equal to 2 

MM (i?) = max(0, min (R, 2)) 


• the van Leer limiter 


■9 vl (R) 


R + \R | 

1 + 1*1 ’ 


Other limiter formulations involving three successive one-dimensional slopes are given in 
Jameson, 1993; Cournede, Debiez and Dervieux, 1998. 


4. Further Advanced Topics 

The material presented in previous sections gives a brief overview of the derivation and analysis 
of finite volume methods. For simplicity and brevity of the presentation, exclusive attention has 
been devoted to scalar nonlinear conservation laws in divergence form. In this overview, special 
consideration has been given to the formulation and stability analysis of higher order accurate 
schemes since these developments have had the largest impact on development of industrial 
computer codes in use at the time of this writing. The remainder of this chapter will consider 
several extensions of the finite volume method. Section 4.1 considers a class of higher order 
accurate discretizations in time that still preserve the stability properties of the fully-discrete 
schemes using Euler time integration. This is followed by a discussion of generalizations of the 
finite volume method for problems including second order diffusion terms and the extension 
to systems of nonlinear conservation laws. 

4.1 . Higher order time integration schemes 

The derivation of finite volume schemes in Sect. 2 first yielded a semi-discrete formulation 
(21) that was later extended to a fully-discrete formulation (22) by the introduction of first 
order accurate forward Euler time integration. These latter schemes where then subsequently 
extended to higher order accuracy in space using a variety of techniques. For many computing 
problems of interest, first order accuracy in time is then no longer enough. To overcome this 
low order accuracy in time, a general class of higher order accurate time integration methods 
was developed that preserve stability properties of the fully-discrete scheme with forward Euler 
time integration. Following Gottlieb, Shu and Tadmor, 2001 and Shu, 2001, these methods will 
be referred to as Strong Stability Preserving (SSP) time integration methods. 

Explicit SSP Runge-Kutta methods were originally developed by Shu, 1988; Shu and Osher, 
1988 and Gottlieb and Shu, 1998 and called TVD Runge-Kutta time discretizations. In a 
slightly more general approach, total variation bounded (TVB) Runge-Kutta methods were 
considered by Cockburn and Shu, 1989; Cockburn, Lin and Shu, 1989; Cockburn, Hou and Shu, 
1990; Cockburn and Shu, 1998 in combination with the discontinuous Galerkin discretization 
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in space. Kiither, 2000 later gave error estimates for second order TVD Runge-Kutta finite 
volume approximations of hyperbolic conservation laws. 

To present the general framework of SSP Runge-Kutta methods, consider writing the semi- 
discrete finite volume method in the following form 

fu(t) = L(U(t)) (93) 

where U — U(t) denotes the solution vector of the semi-discrete finite volume method. Using 
this notation together with forward Euler time integration yields the fully-discrete form 

U n+ 1 = U n -At L(U n ) , (94) 

where U n is now an approximation of U(t n ). As demonstrated in Section 2.2, the forward 
Euler time discretization is stable with respect to the L°°-norm, i.e. 

||*7 n+1 ||oo < imioo . (95) 

subject to a CFL-like time step restriction 

At<At 0 • (96) 

With this assumption, a time integration method is said to be SSP (see Gottlieb, Shu and 
Tadmor, 2001) if it preserves the stability property (95), albeit with perhaps a slightly different 
restriction on the time step 

At < c At 0 (97) 

where c is called the CFL coefficient of the SSP method. In this framework, a general objective 
is to find SSP methods that are higher order accurate, have low computational cost and storage 
requirements, and have preferably a large CFL coefficient. Note that the TVB Runge-Kutta 
methods can be embedded into this class if the following relaxed notion of stability is assumed 

||C/ n+1 || 00 < (1 + O(At)) ||C/ n || 00 . (98) 

4-1-1- Explicit SSP Runge-Kutta methods . Following Shu and Osher, 1988 and the review 
articles by Gottlieb, Shu and Tadmor, 2001; Shu, 2001, a general m stage Runge-Kutta method 
for integrating (93) in time can be algorithmically represented as 

U° := U n , 

u l := A* £(£/*)), a lk > 0, l = (99) 

fc=0 

jyn+1 jjvn 

To ensure consistency, the additional constraint ^i=o a ik = 1 is imposed. If, in addition, all 
fiat are assumed to be non-negative, it is straightforward to see that the method can be written 
as a convex (positive weighted) combination of simple forward Euler steps with At replaced 
by §^At. From this property, Shu and Osher, 1988 concluded the following lemma: 

Lemma 4.1. If the forward Euler method (94) is L°° -stable subject to the CFL condition (96), 
then the Rung-Kutta method (99) with f3ik > 0 is SSP, i.e. the method is L°° -stable under the 
time step restriction (97) with CFL coefficient 

c = min . (100) 

<*lk 

In the case of negative a similar result can be proven, see Shu and Osher, 1988. 
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4.1.2. Optimal second and third order nonlinear SSP Runge-Kutta methods. Gottlieb, Shu 
and Tadmor, 2001 [Proposition 3.1] show that the maximal CFL coefficient for any m-stage, m- 
th order accurate SSP Runge-Kutta methods is c = 1. Therefore, SSP Runge-Kutta methods 
that achieve c = 1 are termed “optimal”. Note that this restriction is not true if the number 
of stages is higher than the order of accuracy, see Shu, 1988. 

Optimal second and third order nonlinear SSP Runge-Kutta methods are given in Shu and 
Osher, 1988. The optimal second order, two-stage non-linear SSP Runge-Kutta method is 
given by 

U° := U n , 

U 1 := U°+AtL(U°) , 

U n+1 := \u Q + \iJ 1 + \At L{U X ) . 

£ £ £ 

This method corresponds to the well known method of Heun. Similarly, the optimal third 
order, three-stage non-lineax SSP Runge-Kutta method is given by 

U° := U n , 

U 1 := U° + At L(U°) , ■ 

U 2 := -V° + -V 1 + HU 1 ) , 

4 4 4 

u n+1 := \ir° + \u 2 + \At L{U 2 ) . 

000 

Further methods addressing even higher order accuracy or lower storage requirements are given 
in the review articles of Gottlieb, Shu and Tadmor, 2001 and Shu, 2001 where SSP multi-step 
methods are also discussed. 

4.2. Discretization of elliptic problems 

Finite volume methods for elliptic boundary value problems have been proposed and analyzed 
under a variety of different names: box methods, covolume methods, diamond cell methods, 
integral finite difference methods and finite volume element methods, see Bank and Rose, 
1987; Cai, 1991; Siili, 1991; Lazarov, Michev and Vassilevsky, 1996; Viozat et al., 1998; 
Chatzipantelidis, 1999; Chou and Li, 2000; Hermeline, 2000; Eymard, Galluoet and Herbin, 
2000; Ewing, Lin and Lin, 2002. These methods address the discretization of the following 
standard elliptic problem in a convex polygonal domain Q, C R 2 

-V- AVu = / in n , ' (101) 

u(x ) = 0 on dQ 

for A G l 2x2 , a symmetric positive definite matrix (assumed constant). Provided / E H&(Q) 
then a solution u exists such that u E i?o +2 (n), — 1 < /? < 1, /? ^ ±1/2, where H s ( fi) denotes 
the Sobolev space of order s in fi. 

Nearly all the above mentioned methods can be recast in Petrov-Galerkin form using a 
piecewise constant test space together with a conforming trial space. A notable exception is 
given in Chatzipantelidis, 1999 wherein nonconforming Crouzeix-Raviart elements are utilized 
and analyzed. To formulate and analyze the Petrov-Galerkin representation, two tessellations 
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of ft are considered: a triangulation T with simplicial elements K € T and a dual tessellation 
T * with control volumes T € T*. In the class of conforming trial space methods such as the 
finite volume element (FVE) method, a globally continuous, piecewise p-th order polynomial 
trial space with zero trace value on the physical domain boundary is constructed 

x h = {v e c°((l) | v| K e V P (K), VK € T and uj an = 0} 

using nodal Lagrange elements on the simplicial mesh. A dual tessellation T* of the Lagrange 
element is then constructed, see for example Fig. 10 which shows a linear Lagrange element 
with two dual tessellation possibilities. These dual tessellated regions form control volumes 
for the finite volume method. The tessellation technique extends to higher order Lagrange 



a. Voronoi dual tessellation b. Median dual tessellation 


Figure 10, Two control volume variants used in the finite volume discretization of second order 
derivative terms: (a) Voronoi dual where edges of the Voronoi dual are perpendicular to edges 
of the triangulation and (b) median dual formed from median dual segments in each triangle. 

elements in a straightforward way. A piecewise constant test space is then constructed using 

r* 

Y h = {v\v\ T e x (T), vrer*} 

where x(T) is a characteristic function in the control volume T. The finite volume element 
discretization of (101) then yields the following Petrov-Galerkin formulation: Find Uh € Xh 
such that 

Y] { [ Wh AVu h -dv+ f w h fdx)= 0 , Vw h £ Y h . (102) 

vrer* ^ JdT Jt ' 

The analysis of (102) by Ewing, Lin and Lin, 2002 using linear elements gives an a priori 
estimate in an L 2 norm that requires the least amount of solution regularity when compared 
to previous methods of analysis. 

Theorem 4.2 (FVE a priori error estimate, Ewing, Lin and Lin, 2002) Assume a 2- 
D quasi-uniform triangulation T with dual tessellation T* such that 3 C > 0 satisfying 

C^h 2 <\T\<Ch 2 , VT€T* . 

Assume that u and Uh are solutions of (101) and (102) respectively with u € H 2 (Q), 
f € H& , (0 < /? < 1). Then 3C f > 0 such that the a priori estimate holds 

ll w ~ Uh\\L 2 (Q)-< C' {h 2 |M|jy2(n) +/l 1+ ^ ||/||jf*(«)) • (103) 
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Unlike the finite element method, the error estimate (103) reveals that optimal order 
convergence is obtained only if / G H& with /3 > 1. Moreover, numerical results show that the 
source term regularity can not be reduced without deteriorating the measured convergence rate. 
Optimal convergence rates are also shown for the nonconforming Crouzeix-Raviart element 
based finite volume method analyzed by Chatzipantelidis, 1999 for u € H 2 (Q) and / G 
An extensive presentation and analysis of finite volume methods for elliptic equations 
without utilizing a Petrov-Galerkin formulation is given in Eymard, Galluoet and Herbin, 2000. 
In this work, general boundary conditions that include non-homogeneous Dirichlet, Neumann 
and Robin conditions are discussed. In addition, the analysis is extended to general elliptic 
problems in divergence form including convection, reaction and singular source terms. 


4.3. Conservation laws including diffusion terms 

As demonstrated in Sect. 1, hyperbolic conservation laws are often approximations to physical 
problems with small or neaxly vanishing viscosity. In other problems, the quantitative solution 
effects of these small viscosity terms are actually sought. Consequently, it is necessary in these 
problems to include viscosity terms into the conservation law formulation. As a model for these 
latter problems, a second order Laplacian term with small diffusion parameter is added to the 
first order Cauchy problem, i.e. 


dtu + V • f(u) - eAu = 0 inE d xR + , (104a) 

u(x, 0) = uq in R d . (104b) 

Here u(x,t) : R d x R+ -► E denotes the dependent solution variable, / G C( M) the hyperbolic 
flux and e > 0 a small diffusion coefficient. Application of the divergence and Gauss theorems 
to (104a) integrated in a region T yields the following integral conservation law form 



• dv — 



eVu • dv = 0 . 


(105) 


A first goal is to extend the fully-discrete form (22) of Sect. 2 to the integral conservation law 
(105) by the introduction of a numerical diffusion flux function djk(uh ) for a control volume 
Tj G T such that 

eVu • dv « £ djk{u h ) . 

e,fc €dTj 

When combined with the general finite volume formulation (22) for hyperbolic conservation 
laws, the following fully-discrete scheme is produced 

A fTl 

u ? +1 = u j - fpj E (9 jk (u?,uT)-d jk m) , VTj e T . (106) 

1 e jk edTj 

In this equation, the index m may be chosen either asnorn+1, corresponding to an explicit 
or implicit discretization. 



4.3.1. Choices of the numerical diffusion flux djk . The particular choice of the numerical 
diffusion flux function djk depends on the type of control volume that is used. Since the 
approximate solution Uh is assumed to be a piecewise constant function, the definition of djk 
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involves a gradient reconstruction of uh in the normal direction to each cell interface ej^. The 
reconstruction using piecewise constant gradients is relatively straightforward if the control 
volumes are vertex-centered, or if the cell interfaces are perpendicular to the straight lines 
connecting the storage locations (see Fig. 10). 

Vertex- centered finite volume schemes. In the case of vertex-centered control volumes such as 
the median dual control volume, a globally continuous, piecewise linear approximate solution 
Uh is first reconstructed on the primal mesh. Vuh is then continuous on the control volume 
interfaces and the numerical diffusion flux straightforwardly computed as 

d jk {v%)= f Vu 1 h-dv jk . (107) 

de jk 


Cell- centered finite volume schemes . In the case of cell-centered finite volume schemes where 
an underlying primal-dual mesh relationship may not exist, a simple numerical diffusion flux 
can be constructed whenever cell interfaces are exactly or approximately perpendicular to the 
straight lines connecting the storage locations, e.g. Voronoi meshes, quadrilateral meshes, etc. 
In these cases, the reconstructed gradient of Uh projected normal to the cell interface ejk can 
be represented by 


Vu™ - v jk 


u™ - uf 
\Xk - Xj\ 


where Xi denotes the storage location of cell T{. The numerical diffusion flux for this case is 
then given by 

d jk {uT)^jJ^- { (uT-uf) . (108) 

Further possible constructions and generalizations are given in Eymard, Gallouet and Herbin, 
2001; Gallouet, Herbin and Vignal, 2000; Herbin and Ohlberger, 2002. 


4.3.2 . Note on stability , convergence and error estimates. Stability analysis reveals a CFL- 
like stability condition for the explicit scheme choice (m = n) in (106) 


A t n < 


a 3 (/Cun ) 2 

aL 9 h min+ £ 


where L g denotes the Lipschitz constant of the hyperbolic numerical flux, a is a positive mesh 
dependent parameter and e is the diffusion coefficient. In constructing this bound, a certain 
form of shape regularity is assumed such that there exists an a > 0 such that for all j, k with 
hk = diam (T*) 

Cih\ < \T k \, a\dT k \<h k , ah k <\x k -xi\. (109) 

Thus, A t n is of the order h 2 for large e and of the order h for e < h. In cases where the 
diffusion coefficient is larger than the mesh size, it is therefore advisable to use an implicit 
scheme (m = n -F 1). In this latter situation, no time step restriction has to be imposed (see 
Eymard et al., 2002; Ohlberger, 2001b). 

In order to demonstrate the main difficulties when analyzing convection dominated problems, 
consider the following result from Feistauer et al., 1999 for a homogeneous diffusive boundary 
value problem. In this work, a mixed finite volume finite element method sharing similarities 
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with the methods described above is used to obtain a numerical approximation Uh of the 
exact solution u. Using typical energy-based techniques, they prove the following a priori 
error bound. 

Theorem 4.3. For initial data u 0 £ L°°(R 2 ) fl PU 1,2 (R 2 ) and r > 0 there exist constants 
Ci , C 2 >0 independent of e such that 

IK->* n ) - «/»(•> * n )IU 2 (n) <cihe C2T/t . (110) 

This estimate is fundamentally different from estimates for the purely hyperbolic problems of 
Sects. 2 and 3. Specifically, this result shows how the estimate strongly depends on the small 
parameter e; ultimately becoming unbounded as e tends to zero. 

In the context of convection dominated or degenerate parabolic equations, Kruzkov- 
techniques have been recently used by Carrillo, 1999; Karlsen and Risebro, 2000 in proving 
uniqueness and stability of solutions. Utilizing these techniques, convergence of finite volume 
schemes (uniform with respect to e -> 0) was proven in Eymard et al., 2002 and a priori 
error estimates were obtained for viscuous approximations in Jakobsen and Karlsen, 2001 and 
Eymard, Gallouet and Herbin, 2002. Finally, in Ohlberger, 2001a; Ohlberger, 2001b uniform 
a posteriori error estimates suitable for adaptive meshing are given. 

4-4- Extension to systems of nonlinear conservation laws 

A positive attribute of finite volume methods is the relative ease in which the numerical 
discretization schemes of Sects. 2 and 3 can be algorithmically extended to systems of nonlinear 
conservation laws of the form 

dtu + V * f{u) = 0 in ]R d x M + , (111a) 

u(x,0) = uo(x) in R d (111b) 

where u(x,t) : R d x RA — > R m denotes the vector of dependent solution variables, f(u) : 
R m i — y R mxd denotes the flux vector, and uo(x) : R d -> R m denotes the initial data vector at 
time t = 0. It is assumed that this system is strictly hyperbolic, i.e. the eigenvalues of the flux 
Jacobian A(v; u) B df /du • v are real and distinct for all bounded v £ R d . 

The main task in extending finite volume methods to systems of nonlinear conservation laws 
is the construction of a suitable numerical flux function. To gain insight into this task, consider 
the one-dimensional linear Cauchy problem for u(x,t) :lxl + 4 R m and uq(x) : R R m 

dtu + d x (Au) = 0 in R x R + , 

u(z,0) = Uq(x) in R (112) 

where A £ R mxm is a constant matrix. Assume the matrix A has m real and distinct 
eigenvalues, Ai < X\ < • • * < A m , with corresponding right and left eigenvectors denoted 
by Tk € R m and Ik 6 M m respectively for k = 1, . . . ,m. Furthermore, let X £ K mxm denote 
the matrix of right eigenvectors, X — [tt, . . . ,r m ], and A £ M mxm the diagonal matrix of 
eigenvalues, A = diag(Ai, . . . , A m ) so that A = XAX~ l . The one-dimensional system (112) 
is readily decoupled into scalar equations via the transformation into characteristic variables 
a = X^u 


d t a + d x (Aa) = 0 in R x R + , 
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a(z,0) — a 0 (x) in 1 


( 113 ) 


and component-wise solved exactly 

a ^ (x, t) = a ^ (x- X k t) , fc = l,...,m 
or recombined in terms of the original variables 


m 

u(x, t) = ^2 l k ■ U 0 (x - Afc t) r k . 

k = 1 


Using this solution, it is straightforward to solve exactly the associated Riemann problem for 
w(£,r) G R 171 

d T w + d^(Aw) = 0 in R x E + 


with initial data 



if e < o 

if ^ > 0 


thereby producing the following Godunov-like numerical flux function 


g(u,v) = Aw(0, R + ) 

= \{Au + Av)- l -\A\ (v-u) 


(114) 


with |A| = X|A|X _1 . When used in one-dimensional discretization together with piecewise 
constant solution representation, the linear numerical flux (114) produces the well-known 
Courant-Isaacson-Rees (CIR) upwind scheme for linear systems of hyperbolic equations 


= «* - £ ( A+ W - 1) + (<•?+. - <)) 

where A ± = XA ± X~ 1 . Note that higher order accurate finite volume methods with slope 
limiting procedures formally extend to this linear system via component wise slope limiting of 
the characteristic components , k = 1, . . . m for use in the numerical flux (114). 


4-4-1 . Numerical flux functions for systems of conservation laws . In Godunov’s original work 
(see Godunov, 1959), exact solutions of the one-dimensional nonlinear Riemann problem of 
gas dynamics were used in the construction of a similar numerical flux function 

g G (u,v) = f(w( 0,R + ))-z/ (115) 


where w(f,r) € W 71 is now a solution of a nonlinear Riemann problem 


d T w A- dzfW(w) = 0 in R x R + 


with initial data 



if e < o 

if £ > 0 ‘ 


Recall that solutions of the Riemann problem for gas dynamic systems are a composition of 
shock, contact and rarefaction wave family solutions. For the gas dynamic equations considered 
by Godunov, a unique solution of the Riemann problem exists for general states u and v 
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except those states producing a vacuum. Even so, the solution of the Riemann problem is 
both mathematically and computationally nontrivial. Consequently, a number of alternative 
numerical fluxes have been proposed that are more computationally efficient. These alternative 
numerical fluxes can be sometimes interpreted as approximate Riemann solvers. A partial list 
of alternative numerical fluxes is given here. A more detailed treatment of this subject is given 
in Godlewski and Raviart, 1991, Kroner, 1997, and LeVeque, 2002. 

• Osher-Solomon flux (Osher and Solomon, 1982). This numerical flux is a system 
generalization of the Enquist- Osher flux of Sect. 2. All wave families are approximated in 
state space as rarefaction or inverted rarefaction waves with Lipschitz continuous partial 
derivatives. The Osher-Solomon numerical flux is of the form 

5 os (u,t;) = ^(/(u) +/(«))• v - ^ J \A{u-w)\dw 

where \A\ denotes the usual matrix absolute value. By integrating on m rarefaction 
wave integral subpaths that are each parallel to a right eigenvector, a system decoupling 
occurs on each subpath integration. Furthermore, for the gas dynamic equations with 
ideal gas law, it is straightforward to construct m — 1 Riemann invariants on each subpath 
thereby eliminating the need for path integration altogether. This reduces the numerical 
flux calculation to purely algebraic computations with special care taken at sonic points, 
see Osher and Solomon, 1982. 

• Roe flux (Roe, 1981). Roe’s numerical flux can be a interpreted as approximating all 
waves families as discontinuities. The numerical flux is of the form 

g Roe (u,v) = i(/(u) + /(«)) -v- ^\A(v;u,v)\(v -u) 

where A( v\ u, v) is the “Roe matrix” satisfying the matrix mean value identity 

(f(v) - /(«)) • V = A{ V] u, v ) ( V - u ) 

with A(i /; ti, u) = A(i /; u). For the equations of gas dynamics with ideal gas law, the Roe 
matrix takes a particularly simple form. Steady discrete mesh-aligned shock profiles are 
resolved with one intermediate point. The Roe flux does not preclude the formation of 
entropy violating expansion shocks unless additional steps are taken near sonic points. 

• Steger- Warming flux vector splitting (Steger and Warming, 1981). Steger and Warming 
considered a splitting of the flux vector for the gas dynamic equations with ideal gas 
law that exploited the fact that the flux vector is homogeneous of degree one in the 
conserved variables. From this homogeneity property, Euler’s identity then yields that 
f(u) • v = A(i/]u) u. Steger and Warming then considered the matrix splitting 

A = A+ + A~ , A ± = XA ± X~ 1 

where is computed component wise. From this matrix splitting, the final upwind 
numerical flux function was constructed as 

g sw (u 7 v) = A* (v\ u) u + A~(i/]v) v . 

Although not part of their explicit construction, for the gas dynamic equations with ideal 
gas law, the jacobian matrix <9<? sw jdu has eigenvalues that are all nonnegative and the 
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Jacobian matrix dg sw jdv has eigenvalues that are all nonpositive whenever the ratio of 
specific heats 7 lies in the interval [1,5/3]. The matrix splitting leads to numerical fluxes 
that do not vary smoothly near sonic and stagnation points. Use of the Steger- Warming 
flux splitting in the schemes of Sect. 2 and 3 results in rather poor resolution of linearly 
degenerate contact waves and velocity slip surfaces due to the introduction of excessive 
artificial diffusion for these wave families. 

• Van Leer flux vector splitting. Van Leer, 1982 provided an alternative flux splitting for 
the gas dynamic equations that produces a numerical flux of the form 

g Vh (u,v) = f~{u) + f + (v) 


using special Mach number polynomials to construct fluxes that remain smooth near 
sonic and stagnation points. As part of the splitting construction, the jacobian matrix 
dg sw /du has eigenvalues that are all nonnegative and the matrix dg sw /dv has 
eigenvalues that are all nonpositive. The resulting expressions for the flux splitting 
are somewhat simpler when compared to the Steger- Warming splitting. The van Leer 
splitting also introduces excessive diffusion in the resolution of linearly degenerate contact 
waves and velocity slip surfaces. 

• System Lax-Friedrichs flux. This numerical flux is the system equation counterpart of 
the scalar Lax-Friedrichs flux (27). For systems of conservation laws the Lax-Friedrichs 
flux is given by 

g LF (u, v ) = !(/(«) + f(v)) ■ v - lc*( v) ( v - u ) 
where a(v) is given through the eigenvalues A* of A(i/]w) 


<*(") = max sup |A*(i/;w)| . 

l < k <™ we[u,v] 


The system Lax-Friedrichs flux is usually not applied on the boundary of domains since it 
generally requires an over specification of boundary data. The system Lax-Friedrichs flux 
introduces a relatively large amount of artificial diffusion when used in the schemes of 
Sect. 2. Consequently, this numerical flux is typically only used together with relatively 
high order reconstruction schemes where the detrimental effects of excessive artificial 
diffusion are mitigated. 

• Harten-Lax-van Leer flux (Harten, Lax and van Leer, 1983). The Harten-Lax-van Leer 
numerical flux originates from a simplified two wave model of more general m wave 
systems such that waves associated with the smallest and largest characteristic speeds 
of the m wave system are always accurately represented in the two wave model. The 
following numerical flux results from this simplified two wave model 


g HLL (u,v) = !(/(«) + f(v))-v 


2 o: max oc m [ n 


Cf(t>) -/(«))• * + 


^max £*min 


(v ~ u) 


where 


ttmaxM = max (0, Sup Ajfe(i/;tu)) , a min = min (0, inf A k (u\w)) . 

l<fc<m w e[u,v] 1 <k<m we[u y v] 

When compared to the Lax-Friedrichs flux, this flux can be considerably more accurate 
in flow situations where 0 < |(a max + a m i n )/(a m ax “ a min)| < 1. Using this flux, full 
upwinding is obtained for supersonic flow. Modifications of this flux are suggested in 
Einfeldt et ah, 1998 to improve the resolution of intermediate waves as well. 
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Further examples of numerical fluxes (among others) include the kinetic flux vector splitting 
due to Deshpande, 1986, the advection upstream splitting flux (AUSM) of Liou and Steffen, 
1993, and the convective upwind and split pressure (CUSP) flux of Jameson, 1993 and Tatsumi 
et al., 1994. 


5. Concluding Remarks 

The literature associated with the foundation and analysis of the finite volume methods is 
extensive. This article gives a very brief overview of finite volume methods with particular 
emphasis on theoretical results that have significantly impacted the design of finite volume 
methods in everyday use at the time of this writing. More extensive presentations and 
references on various topics in this article can be found in the books by Godlewski and Raviart, 
1991, Kroner, 1997, Eymard, Galluoet and Herbin, 2000 and LeVeque, 2002. 
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